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A first  principles  description  of  electronic  excitation  and  spin-orbit  recoupling  in 
alkali-rare-gas  atom  colliding  pairs  is  developed  introducing  /-dependent  pseudopo- 
tentials and  including  two  and  three-body  polarization  terms  and  the  spin-orbit 
interatomic  potential.  The  treatment  combines  an  eikonal  (or  short  wavelength) 
approximation  for  the  nuclear  motion  and  time-dependent  molecular  orbitals  to 
provide  interatomic  potentials,  their  non-adiabatic  couplings,  and  state  populations 
during  interactions. 

Our  pseudopotential  choice  is  adequate  for  these  systems  in  terms  of  accuracy.  A 
description  of  the  matrix  equations  and  computational  details  of  the  pseudopotential 
code  are  presented.  Results  on  ionization  energies  of  the  Li  and  Na  atoms  and  the 
potential  energy  functions  of  distance  for  LiHe,  LiNe,  NaHe  and  NaNe  are  presented. 
Our  results  are  in  excellent  agreement  with  experiment  and  other  theories. 

The  theory  of  the  eikonal  time-dependent  molecular  orbital  (Eik/TDMO)  method 
in  terms  of  pseudopotentials  is  presented  as  well  as  a detailed  description  of  the  com- 
putational approach.  We  study  the  time-dependence  of  orbital  alignment,  orienta- 
tion and  population.  We  discuss  the  effects  of  the  basis  set  size  on  the  calculations, 
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and  compare  our  results  with  experiment  and  other  calculations.  Our  integral  cross 
sections  for  LiHe,  LiNe,  NaHe  and  NaNe,  obtained  with  a large  basis  set,  are  in 
excellent  agreement  with  experiment. 

The  theory  for  recoupling  of  angular  momenta  in  alkali-rare-gas  atom  thermal 
collisions  is  developed  and  the  computational  aspects  of  the  spin-orbit  coupling  code 
in  terms  of  pseudopotentials  are  presented.  Results  for  the  spin-orbit  splitting  of 
the  2P  states  of  Li  and  Na  are  presented  along  with  fine-structure  cross  sections  for 
collisions  at  400  and  450  K.  The  agreement  with  experiment  and  other  theories  is 
very  good. 
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CHAPTER  1 
INTRODUCTION 

This  study  is  part  of  a broader  effort  to  bring  new  insights  into  the  time- 
dependence  of  collision  of  quasi-one-electron  systems.  We  are  particularly  interested 
in  studying  the  time-dependence  of  electronic  excitation  and  the  time-dependence 
of  spin-orbit  coupling  in  slow  alkali-rare  gas  atom  collisions. 

The  work  developed  in  this  dissertation  is  related  to  phenomena  of  experimental 
interest,  such  as  pressure  broadening  of  spectra  in  condensed  matter  and  kinetics 
in  hot  gases  and  plasmas.1-5  In  a related  area,  the  optical  spectroscopy  of  alkali 
atoms  in  rare  gas  atom  clusters  has  been  the  focus  of  recent  research,  particularly 
in  connection  with  spectra  in  liquid  helium.6-9 

It  is  of  conceptual  value  to  gain  a fundamental  understanding  of  the  most  basic 
time-dependent  aspects  of  the  rearrangement  phenomena  in  simple  systems  where 
detailed  studies  can  be  carried  out. 

Alkali-rare  gas  atom  collisions  constitute  a particularly  important  subject  of 
study.  In  addition  to  the  fundamental  interesting  questions  they  pose,  they  are  at 
the  heart  of  a large  number  of  applications  in  such  areas  as  laser  systems,  surface 
interactions,  chemical  reactions,  energy  transfer  and  others.10 

In  particular,  we  are  interested  here  in  studying  the  time-evolution  of  electronic 
excitation  and  spin-orbit  recoupling  of  Li  and  Na  in  collisions  with  He  and  Ne  in 
terms  of  a formalism  developed  to  consider  our  many  electron  system  as  a one-active 
electron  system. 

1.1  Overview  of  Experimental  Work 

Over  the  last  thirty  years,  several  experimental  results  for  these  systems  have 
improved  our  current  understanding  of  the  electronic  rearrangement  phenomena. 
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Overviews  of  those  studies  are  available  in  the  literature.11,12  Experimentalists  have 
been  able  to  select  the  initial  states  of  colliding  species  to  perform  collisions  at  spe- 
cific energies.  The  level  of  control  in  experiments  has  reached  the  point  at  which 
theoretical  results  from  first  principles  calculations  can  be  of  help  in  the  interpreta- 
tion of  experiments.5 

Excitation  of  the  alkali  projectile  in  the  energy  range  from  a fraction  of  a keV 
to  100  keV  has  been  investigated.  Representative  studies  for  emission  cross  sections, 
alignment  and  polarization13-16  have  partly  revealed  the  mechanisms  responsible  for 
the  observed  excitations. 

The  literature  on  experimental  aspects  of  spin-orbit  effects  in  gas-phase  chem- 
ical reactions  done  before  mid-80’s  has  been  reviewed.17'18  More  recently,  some 
experimental  work  has  been  done  to  understand  the  role  of  spin-orbit  coupling  in 
the  stereochemical  behavior  of  larger  systems.  The  use  of  spin-polarized  beams  is 
bringing  new  insights  because  collisions  can  be  studied  in  specific  spin  states.  The 
use  of  coincidence  detection  allows  observation  of  a larger  number  of  parameters. 
The  ability  to  use  the  polarization  of  lasers  to  prepare  highly  aligned,  or  direc- 
tional, reagents  has  opened  new  possibilities  in  the  area  of  orbital  stereochemistry. 
Recent  advances  in  laser  pump-probe  techniques  have  been  applied  to  the  study 
fine-structure  changing  in  collisions  of  alkali  atoms  with  rare  gases.19,20  Important 
advances  in  experimental  spectroscopic  studies  of  alkali  atoms  attached  to  large  rare- 
gas  atom  clusters  have  allowed  probing  and  modeling  as  the  alkali  atom  is  excited 
or  deexcited.6-9 

1.2  Overview  of  Theory 

The  field  of  quantum  molecular  dynamics  has  seen  unparalled  development  over 
the  last  three  decades.  Methods  for  the  solution  of  the  time-dependent  Schrodinger 
equation  for  atomic  collisions  process  abound  in  the  literature.  However  the  ma- 
jority of  these  methods  make  use  of  the  central  dogma  in  both  time-dependent  and 
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time-independent  calculation  of  molecular  properties  and  reaction  rates:  the  Born- 
Oppenheimer  (BO)  approximation.  The  BO  approximation  is  invoked  to  separate 
the  Schrodinger  equation  for  the  whole  system  into  two  separate  equations  to  be 
solved  independently:  one  for  the  electrons  and  another  one  for  the  nuclei.  Time- 
dependent  methods  can  be  classified  in  two  groups  based  on  the  way  that  they 
approximate  the  molecular  time-dependent  Schrodinger  equation.  The  first  group 
consists  of  methods  that  solve  the  time-dependent  Schrodinger  equation  for  the  nu- 
clear degrees  of  freedom  subject  to  known  potential  energy  surfaces  (PES),  while 
in  the  second  group  the  electronic  and  nuclear  degrees  of  freedom  are  propagated 
simultaneously.  This  is  different  from  finding  the  solution  for  the  electronic  prob- 
lem along  a nuclear  trajectory.  The  eikonal  time-dependent  Hartree-Fock  method 
(Eik/TDHF)  method21-24  and  the  electron  nuclear  dynamics  method  (END)25  can 
be  cataloged  in  the  second  group. 

Because  of  the  success  of  methods  in  the  first  group  for  the  study  of  collisional 
systems,  explicit  time-dependent  calculations  for  both  electrons  and  nuclei  have  not 
become  common  practice.  However,  there  is  a lot  to  gain  from  these  approaches  as 
a much  clearer  picture  of  excitation  mechanism  can  be  drawn  from  time-dependent 
simulations.  For  example,  in  collisions  involving  excitation  or  charge  transfer,  fol- 
lowing how  charge  changes  with  nuclear  configuration  over  time  can  offer  insight 
into  mechanisms. 

First  principles  methods  of  full  dynamics  have  some  obvious  advantages  when 
compared  to  other  methods:  there  is  no  need  to  use  inconvenient  coordinate  systems 
because  atomic  cartesian  coordinates  are  used  throughout  and  there  is  also  no  need 
to  construct  potential  energy  surfaces. 

1.3  Overview  of  Theoretical  Studies  in  Alkali- Rare-Gas  Interactions 

A complete  and  rigorous  theory  of  these  processes  involves  calculation  of  the  full 
quantum  wave  function  ip  describing  the  motion  of  the  nuclei  and  active  electron 
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with  solutions  of  close-coupling  scattering  equations.  This  formulation  has  been 
implemented  by  different  authors.4, 26-31 

The  atomic-state  expansion  method  was  used  by  Andersen  and  Nielsen32  for 
Li  and  Na  with  He  and  Ne  collision.  There,  they  adopted  various  semiempirical 
model  potentials  or  the  pseudopotential  proposed  by  Bottcher33  and  Baylis34  to  ob- 
tain electronic  wavefunction  with  which  to  carry  out  the  close-coupling  calculations. 
They  also  used  a frozen  core  Hartree-Fock  potential.  Their  calculated  integral  cross 
sections  agree  well  with  experiment.  However  the  shape  of  the  cross  sections  differs 
from  the  experimental  measurements  in  the  low-energy  region. 

More  recent  calculations  by  Kimura  and  Pascale35  present  results  for  the  integral 
cross  sections  of  the  alkali  atom-He  system  in  the  low  keV  energy  regime.  Their 
results  are  obtained  with  a close-coupling  method  including  electron  translation 
factors  and  straight  line  trajectories.  The  agreement  with  experiment  is  good  for 
LiHe  and  NaHe.  For  KHe  and  CsHe  the  results  are  only  qualitatively  good.  More 
recently  Archer  et  al2  presented  results  on  the  orientation  and  alignment  of  LiHe 
and  NaHe  collisions. 

The  close-coupling  method  is  known  to  give  accurate  cross  sections;  however  it 
has  some  disadvantages  because  the  formulation  of  the  theory  and  the  long  numerical 
codes  are  complex,  and  physical  insight  can  be  lost  when  many  states  are  coupled.  A 
more  intuitive,  though  less  rigorous,  semiclassical  model  has  successfully  been  used 
to  study  these  systems. 

A strong  coupling  approximation  to  obtain  analytic  formulas  for  the  cross  sec- 
tions of  fine-structure  transitions  in  the  Na-Ar  system  is  presented  in  Nikitin.36 
Some  authors  preferred  to  numerically  solve  equations  for  the  Na-He  system  using 
a semiclassical  impact  parameter  method,37,38  later  refined39  to  take  into  account 
trajectory  effects.  More  recently40  a semiclassical  formalism  was  used  to  monitor 
the  expectation  values  of  the  orbital  and  spin  angular  momentum  vectors. 
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1.4  Our  Approach 

The  Eik/TDMO  method  developed  in  our  group  is  a semi-classical  method  that 
uses  a combination  of  the  eikonal  and  the  TDMO  approximations.  Starting  with 
an  eikonal  representation  of  the  total  wave  function,  a wave  function  is  constructed 
from  the  classical  trajectories  and  the  TDMO  formulation  is  developed  in  terms  of  a 
density  operator.  The  coupling  of  the  density  operator  to  nuclear  motions  is  solved 
with  a “relax  and  drive”  procedure21,22  that  couples  the  fast  (electronic)  and  slow 
(nuclear)  degrees  of  freedom.  We  have  made  this  one-electron  method  capable  of 
studying  the  quasi-one  electron  system  where  the  valence  electron  has  been  described 
explicitly  while  its  interaction  with  the  alkali  core  and  the  rare  gas  atom  has  been 
described  with  /-dependent  atomic  pseudopotentials.41,42 

1.5  Systems  To  Be  Studied 

This  dissertation  studies  the  first-principles  molecular  dynamics  of  the  alkali- 
rare-gas  atom  colliding  pair.  We  will  concentrate  on  two  different  energy  regimes: 

A high  energy  regime  (with  collision  energies  between  500  eV  and  30  keV)  in 
which  we  study  the  electronic  excitation  of  the  alkali  atom  M from  the  ground  state 
to  any  excited  state 


M(nl)  + Rg  + Rg, 

where  n and  l are  the  principal  and  orbital  quantum  numbers  and  M = Li,  Na 
and  Rg  = He,  Ne.  We  will  monitor  the  time-dependence  of  the  atomic  populations 
and  orbital  alignment  and  orientation  and  we  will  also  calculate  final  state-to  state 
integral  cross  sections.  The  terms  orientation  and  alignment  refer  to  parameters 
related  to  the  shape  and  dynamics  of  an  excited  atomic  or  molecular  level,  as  it  is 
manifested  in  a non-statistical  population  of  the  magnetic  sublevels. 

Hyperthermal  energy  regime  (with  energies  of  only  a fraction  of  an  eV)  where 
we  will  study  the  recoupling  of  angular  momenta  in  collisions  of  excited  Na  with  He. 
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These  collisions  may  result  either  in  reorientation 

Na*(2PjtMj)  + He^  Na\2PJM .)  + He  (1.1) 

or  in  intramultiplet  transitions 

Na*(2PJMj ) + He^  N a*  {2  Pj,  ^ + He  (1.2) 

with  energy  transfer 

AJ51  = Ej'Mj  — (1.3) 

where  we  have  used  the  labels  J for  the  total  electronic  angular  momentum  and  Mj 
for  its  projection.  At  these  low  energies  these  are  the  only  processes  which  occur 
with  appreciable  probability.11  Spin-orbit  coupling  plays  an  important  role  in  both 
processes.  To  get  a complete  picture  of  the  collisional  processes,  we  we  will  monitor 
the  expectation  values  of  the  electronic  orbital  (1),  spin  (s)  and  total  (j)  angular 
momentum  vectors  along  a trajectory. 

1.6  Outline  of  the  Dissertation 

Chapter  2 gives  a general  introduction  to  the  pseudopotential  theory  and  dis- 
cusses in  more  detail  the  pseudopotential  method  used  in  this  dissertation.  In  that 
chapter,  we  concentrate  on  the  evaluation  of  the  hamiltonian  matrix  elements.  We 
also  discuss  the  effects  of  the  basis  set  size  on  the  calculated  energies.  We  explain 
the  computational  aspects  of  the  pseudopotential  code  and  finally  compare  our  cal- 
culated ionization  energies  for  Li  and  Na  with  other  theories  and  experiments  as 
well  as  the  potential  energy  curves  of  the  alkali-rare-gas  atom  systems:  LiHe,  LiNe, 
NaHe  and  NaNe. 

Chapter  3 generalizes  the  eikonal  time-dependent  molecular  orbital  (Eik/TDMO) 
method23,24  to  include  pseudopotentials  (Eik/TDMO-PP).  Theory  and  computa- 
tional details  are  introduced  in  this  chapter.  The  calculation  of  some  time-dependent 
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physical  properties  such  as  electronic  populations  and  orientation  and  alignment  pa- 
rameters are  reviewed  here  too.  Calculations  of  time-dependent  properties  and  also 
of  the  excitation  integral  cross  sections  are  presented  along  with  their  comparison 
to  experiments  and  other  theories. 

In  Chapter  4,  the  problem  of  recoupling  of  angular  momenta  in  thermal  alkali- 
rare-gas  atom  collisions  is  formulated  within  a formalism  based  on  the  Eik/TDMO- 
PP  method  that  includes  spin-orbit  coupling  (Eik/TDMO-PP-SO).  The  evaluation 
of  the  spin-orbit  coupling  matrix  elements,  and  related  computational  aspects  are 
exposed.  Results  for  the  spin-orbit  coupling  splitting  of  alkali  atoms  and  the  po- 
tential energy  curve  of  the  Na-He  system  are  presented  as  well  as  some  encouraging 
results  on  the  fine-structure  cross  sections  at  400  and  450  K.  Pictures  of  the  time 
evolution  of  the  angular  momenta  vectors  are  presented  here  too. 

Chapter  5 summarizes  the  main  conclusions  obtained  in  this  dissertation. 

Appendix  A presents  the  flow  diagram  of  the  Eik/TDMO  code,  which  includes 
pseudopotentials  and  spin-orbit  coupling.  It  also  contains  some  information  about 
the  main  modules  of  the  program. 

Appendix  B explains  in  detail  the  procedure  we  used  to  solve  the  polarization 
interaction  integrals  of  our  pseudopotential. 

Appendix  C gives  a detailed  explanation  of  the  solution  of  the  short-range 
interaction  integrals  in  terms  of  Cartesian  Gaussian  functions  present  in  our  choice 
of  pseudopotential. 


CHAPTER  2 

PSEUDOPOTENTIAL  THEORY  FOR  ALKALI-RARE-GAS  ATOM  SYSTEMS 

2.1  Introduction 

Much  experimental  and  theoretical  work  has  been  devoted  in  the  last  forty 
years  to  the  study  of  various  processes  occurring  in  collisions  between  ground  state 
or  excited  state  alkali  atoms  and  ground  state  rare-gas  (Rg)  atoms  at  a wide  range 
of  energies.  From  the  experimental  perspective,  the  processes  of  interest  have  been 
inter-doublet  and  intra-doublet  transitions,  quenching,  redistribution,  pressure  spec- 
tral line  broadening,  etc. 

Because  the  interatomic  interactions  are  the  main  physical  quantities  needed  for 
a good  understanding  of  these  collisional  processes,  much  effort  has  been  devoted  to 
calculating  or  experimentally  determining  the  adiabatic  potentials  of  the  alkali-Rg 
systems. 

Standard  ab  initio  all-electron  calculations  of  the  adiabatic  potentials,  which  in 
principle  can  be  very  accurate  when  enough  electronic  configurations  are  included  in 
the  calculation,  become  rapidly  difficult  to  perform  as  the  number  of  electrons  in  the 
atomic  cores  increase.  Calculations  using  pseudopotentials  offer  a very  interesting 
alternative  to  treat  the  problem.  They  take  advantage  of  the  fact  that  the  alkali 
core  and  Rg  atom  have  closed-shell  structures  to  reduce  the  problem  of  the  alkali-Rg 
interaction  to  a three-body  problem. 

This  chapter  gives  a general  introduction  to  the  pseudopotential  theory  and 
discusses  in  detail  the  pseudopotential  method  to  be  used  in  this  work.  It  also 
reports  some  encouraging  results  for  the  ionization  energies  of  Li  and  Na  and  for  the 
potential  energy  curves  (PECs)  of  some  alkali-rare  gas  (alkali-Rg)  diatomic  systems. 
Details  on  the  implementation  of  the  computer  code  are  also  included. 
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2.2  Pseudopotential  Theory 

In  many  problems  of  atomic  and  molecular  physics  the  electrons  of  the  system 
can  be  divided  into  valence  and  core  electrons.  It  is  natural  to  try  to  reduce  those 
problems  from  a mathematical  iV-electron  problem,  where  the  N is  the  total  number 
of  electrons,  to  an  n-electron  problem,  where  n is  the  number  of  valence  electrons. 
Since  in  many  cases  n <C  N,  such  reductions  mean  a significant  conceptual  and 
mathematical  simplification. 

Hellmann,43  Hellmann  and  Kassatotschkin44  and  Gombas45  were  the  first  to 
demonstrate,  using  the  Thomas  -Fermi  model,  that  the  Pauli  exclusion  principle 
for  the  valence  electrons  can  be  replaced  with  an  effective  potential.  Since  the  re- 
quirement that  the  valence  orbital  be  orthogonal  to  the  core  orbital  is  equivalent 
to  the  Pauli  exclusion  principle,  they  replaced  the  need  for  orthogonality  by  the 
effective  potential.  The  introduction  of  the  effective  potential  represents  a con- 
siderable mathematical  simplification,  since  the  orthogonality  requirement  between 
valence  and  core  wave  functions  may  lead  to  difficulties,  especially  in  the  case  of 
large  cores.46 

Later,  in  1959,  Phillips  and  Kleinman47  provided  a rigorous  basis  for  replacing 
the  explicit  core-valence  orthogonality  constraints  by  a modification  of  the  valence 
Hamiltonian.  The  usefulness  of  the  Hellmann  and  the  Phillips-Kleinman  approaches 
is  that  they  provide  the  theoretical  basis  for  the  subsequent  development  of  effective 
potentials. 

From  a mathematical  point  of  view  the  pseudopotential  methods  can  be  divided 
into  two  groups:  ab  initio  and  model  methods.  The  ab  initio  methods  differ  from  the 
model  methods  in  that  parametric  functional  forms  are  not  assumed  for  the  poten- 
tial. Rather,  its  form  is  abstracted  directly  from  the  atomic  Hartree-Fock  equations 
upon  use  of  a specific  pseudo-orbital  transformation.  Examples  for  ab  initio  meth- 
ods are  those  proposed  by  Kahn  et  ah, 48  later  improved  by  Hay  et  al.49  Model 
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pseudopotentials  are  effective  potentials  that  appear  from  simplifications  made  on 
ab  initio  pseudopotentials.  In  general,  they  can  be  semi-empirical  when  they  con- 
tain adjustable  parameters  determined  from  experimental  data,  or  non-empirical 
when  they  do  not  require  experimental  data  for  their  construction,  as  in  the  case  of 
density-dependent  pseudopotentials.50 

Regardless  of  their  differences,  all  the  pseudopotential  methods  have  in  com- 
mon: 

1)  The  potential  should  be  asymptotically  correct  for  large  electron- nucleus  dis- 
tance r,  at  least  to  the  order  of  1/r  and  preferably  higher  orders. 

2)  The  Hamiltonian  operator  eigenvalues  should  accurately  reproduce  the  observed 
spectra  of  the  appropriate  single-valence  electron  system.  There  should  be  no 
eigenvalues  corresponding  to  the  core  orbitals  of  the  physical  system. 

3)  The  potential  should  be  simple  so  that  it  can  be  easily  constructed.  An  analytic 
form  involving  a small  number  of  parameters  is  most  convenient.  This  analytic 
form  should  be  applicable  to  many  atomic  systems,  through  suitable  choice  of 
the  parameters. 

4)  The  potential  should  not  be  large  in  magnitude  in  the  core  region.  This  is 
desirable  to  avoid  computational  difficulties,  and  is  particularly  important  for 
molecular  applications.  In  the  case  of  /-dependent  pseudopotentials,  where  / 
is  the  orbital  angular  momentum  of  the  core  electrons,  for  each  value  of  / one 
introduces  one  or  two  adjustable  parameters  and  so  can  reproduce  exactly  the 
experimental  energy  for  two  levels.  With  two  parameters  for  each  /-value  it  is 
usually  possible  to  obtain  an  excellent  fit  to  all  valence  state  energies. 

The  true  utility  of  effective  potentials  is  that  they  provide  a method  for  reduc- 
ing substantially  the  number  of  calculations,  making  possible  calculation  on  large 
molecules  and  on  molecules  containing  heavy  atoms.  We  have  taken  advantage  of 
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that  feature  in  the  study  of  electronic  excitations  and  electronic  angular  recoupling 
in  collisions  of  alkali  (Li  and  Na)  and  Rg  (He  and  Ne)  atoms. 

Extensive  pseudopotential  molecular-structure  calculations  for  all  the  alkali- 
Rg  systems  have  been  done  using  /-independent  Gombas-type  statistical  pseudopo- 
tentials.34’51-53 Although  these  calculations  yielded  reasonable  potentials  for  the 
systems  involving  heavy  Rg  atoms,  they  fail  badly  in  the  case  of  He  and  Ne.  An  im- 
provement was  achieved  by  applying  the  /-dependent  pseudopotential  technique42,54 
for  LiHe,3,42  LiNe31,41  and  NaNe. 41,54,55  Encouraged  by  the  accuracy  of  those  re- 
sults we  applied  the  /-dependent  pseudopotential  method  for  the  calculation  of  the 
adiabatic  potentials  and  couplings  of  alkali-Rg  systems. 

2.3  Pseudopotential  for  Alkali-Rare  Gas  Atom  Systems 

The  fact  that  an  alkali  core  and  an  Rg  atom  have  closed-shell  structures  allows 
us  to  reduce  a many-body  problem  to  a three-body  problem  (alkali  core,  valence 
electron  and  Rg  atom).  The  interaction  between  the  valence  electron  and  both  the 
alkali  core  and  the  Rg  atom  is  then  represented  by  a pseudopotential  and  if  desired, 
it  can  also  include  polarization  terms  in  order  to  account  for  the  valence  electron- 
core  correlation  effects.  In  the  approach  that  we  have  chosen,  the  effect  of  the  core 
electrons  of  both  the  alkali  atom  and  the  Rg  atom  is  simulated  by  an  /-dependent 
semi-empirical  pseudopotential  including  core  polarization.41  This  pseudopotential 
is  composed  of  two  parts:  the  /-dependent  short-range  term,  and  the  polarization 
potential  containing  a cutoff  parameter. 

2.3.1  Hamiltonian 

Let  us  recall  that  the  alkali-Rg  atomic  pair  is  represented  by  a three-body 
system  (alkali  core,  Rg  atom  and  alkali  valence  electron).  The  Hamiltonian  for 
the  one  electron  alkali-Rg  system  at  a given  internuclear  distance  R and  without 
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spin-orbit  coupling  is  (in  atomic  units) 

HSP  = + vr,  (2.1) 

v%p  = Va(ta)  + Vab(ta,  R.),  (2.2) 

where,  in  general,  r*  ( X — A,  B)  is  the  position  vector  of  the  alkali  valence  electron 
e~  with  respect  to  the  alkali  core  (A)  or  the  Rg  atom  ( B ),  and  R is  the  position 
vector  of  B with  respect  to  A.  The  term  Va{ta)  describes  the  interaction  between 
the  valence  electron  and  the  alkali  core  and  is  given  by 

Va(ta)  = + Vf(TA)  + V}'(rA).  (2.3) 

•A 

Here  ZA  is  the  net  charge  of  the  alkali  core  seen  by  e~  at  an  infinite  distance  and 
Vf  is  the  polarization  potential, 

Vl°\rA)  = ~^wi.rA,  Sc)2,  (2.4) 

where  acA  is  the  dipole  polarizability  of  the  alkali  core  and  w stands  for  a cutoff 
function  defined  by 

w(r,  Sc)  = 1 - exp(— <5cr2),  (2.5) 

with  an  adjustable  cutoff  parameter  Sc.  The  last  term  in  Eq.  2.3  is  the  so-called 
semi-local  (^-dependent)  pseudopotential,  which  is  given  here  as 

VJ>x)  = J>MeXp  (-0,A)Pi.x,  (2.6) 

where  Bu  and  are  pseudopotential  parameters  adjusted  to  experimental  data 
and  Piyx  stands  for  the  projection  operator  on  angular  symmetry  l, 

Pi,x  = ^2  I ^m(rx))(Yim(rx)  | . 


m 


(2.7) 
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The  vector  | Fjm(rx))  represents  a spherical-harmonic  function  centered  on  the 
core  X ( rx  = ^x/fx)  • In  turn,  the  term  VAB  in  Eq.  2.2  describes  the  interaction 
between  the  alkali  atom  and  the  Rg  atom  and  is  represented  as 


VAB(rA,  R)  = VB(rB)  + Vec(rA>  R)  + VCC(R).  (2.8) 


The  term  VB{v b)  (*b  = rA  - R)  describes  the  interaction  between  the  valence 
electron  and  the  Rg  atom  and  is  given  by 

Vb( tb)  = V?{tb)  - || w(rB,  S)4  - ^ W(rB , 5 )6.  (2.9) 

The  parameter  aB',  is  defined  as  aqB'  — aB  — 6/?i,  where  /3i  is  the  dynamical 
correction  to  the  static  quadrupole  polarizability  of  the  Rg  atom.  The  short-range 
term  Vb{tb)  is  defined  as  by  Eq.  2.6.  In  turn,  the  term  Vec  in  Eq.  2.8  is  the 
so-called  cross-term,  which  arises  from  the  polarization  of  B by  both  e~  and  A,  and 
is  given  by 


= + (2.10) 

B B 

where  Pi,  and  P2,  are  the  Legendre  polynomials  and  9 is  the  angle  between  the 
vectors  R and  t#. 

Finally,  the  core-core  interaction  is  assumed  to  have  the  form 


VCC(R)  = V£(R)  - - 


a 


B 


a 


q ' 

B 


2 (R2  + d%)2  2 (P2  + 4)3’ 


(2.11) 


where  the  short-range  part  (R)  can  be  calculated  with  different  methods.34, 42,56,57 
Our  choice  for  the  was  based  on  a more  practical  approach  discussed  later  in 
this  chapter.  The  parameters  aB  , aB'  and  dB  for  He  and  Ne  are  available  in  the 
literature.  In  more  sophisticated  calculations  the  Vcc  term  can  include  exchange 
effects,  Coulomb  interaction  as  well  as  the  dispersion  forces. 


14 


2.3.2  Choice  of  Parameters 

The  pseudopotential  parameters  for  alkali  atoms  employed  for  (Li,  Na)  were 
taken  from  Fuentealba  et  al.58  (see  Table  2-1).  These  parameters  were  determined 
by  adjusting  the  calculated  atomic  energies  to  the  corresponding  experimental  ion- 
ization energies.59  The  pseudopotential  parameters  of  the  Rg  atoms  (He,  Ne)  were 
obtained  from  Czuchaj  et  al.41  (see  Table  2-2).  They  were  determined  by  fitting 
the  experimental  /-wave  scattering  phase  shifts  together  with  the  scattering  length, 
considered  currently  to  be  the  best  available,  for  e_-He  (Ne,  Ar)  elastic  scattering 
in  the  very  low  energy  region. 

The  long-range  polarization  potentials  contain  the  cutoff  functions.  They  are 
simply  expressed  as  suitable  powers  of  the  function  defined  by  Eq.  2.5  with  the 
cutoff  parameter  5 assumed  to  be  the  same  for  both  the  dipole  and  quadrupole 
interaction  terms. 

The  Gaussian  form  of  the  short-range  part  V^r{rB)  of  the  effective  potential 
is  one  of  the  most  widely  encountered  in  the  literature.  Because  the  term  Vg(rB) 
represents  an  effective  potential  acting  between  the  scattered  electron  and  the  atom 
in  the  core  region,  its  precise  form  is  not  essential  for  reproducing  experimental 
scattering  phase  shifts.  Rather  the  form  of  VBr(rB)  is  motivated  by  subsequent  use 
of  the  pseudopotentials  in  molecular  calculations. 

2.3.3  Choice  of  Basis  Sets 

We  have  chosen  Cartesian  Gaussian  functions  <ft*(r;  ctji,k),  j = 1, 2, . . . , JVj, 
to  represent  static  atomic  functions  with  quantum  number  (/)  and  s-,  p-,  and  d- 
symmetry  given  by  a A;  index, 


(2.12) 


k= 1 


providing  segmented  contracted  combinations  of  Gaussians  with  exponents  a and 
coefficients  d obtained  from  the  literature.  There  are  several  choices  of  Cartesian 
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Table  2-1:  Pseudopotential  parameters  for  alkali  atoms  (au) 


aA 

l 

A 

A 

Li 

0.1915 

0.831 

0 

5.786 

1.276 

1 

- 1.065 

1.607 

Na 

0.9947 

0.62 

0 

10.839 

1.378 

1 

2.303 

0.6639 

2 

- 1.777 

0.9249 

K 

5.354 

0.29 

0 

13.564 

0.853 

1 

2.648 

0.3696 

2 

-4.517 

0.6639 

Table  2-2:  Pseudopotential  parameters  for  He  and  Ne  (au) 


nd 

aB 

<4 

A 

<5 

l 

i 

B, 

A 

He 

1.3834 

2.3265 

0.706 

0.79 

0 

1 

0.83 

1.30 

0 

2 

2.27 

0.50 

1 

1 

-0.12 

0.75 

1 

2 

-1.87 

1.00 

Ne 

2.663 

6.458 

1.27 

1.00 

0 

1 

2.50 

0.71 

1 

1 

10.00 

5.20 

1 

2 

1.22 

0.41 

2 

1 

-4.50 

1.00 
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Gaussian  basis  sets  available  in  the  literature,  but  they  are  not  accurate  in  all  of  the 
cases.  The  basis  sets  we  have  used  for  Li  and  Na  were  optimized  in  the  pseudopo- 
tential calculation  of  the  ionization  energies  of  the  ground  and  several  excited  states 
of  s-,  p-  and  d-  symmetries  of  those  alkali  atoms.  For  Li,  we  found  two  finite  basis 
sets,  in  the  usual  notation  (uncontracted/contracted)  Cartesian  Gaussian  functions 
of  s-,  p-  and  d-type,  optimized  for  our  pseudopotential:  a small  (4s4p/2s2p)  basis 
was  taken  from  the  MOLPRO  webpage  and  a big  (9s9p5d/7s7p5d)  basis  taken  from 
Czuchaj  et  al.3  For  Na  we  found  a (13sl2p4d/10s8p4d)  basis  in  Czuchaj  et  al.41 
The  small  basis  sets  obtained  from  the  MOLPRO  webpage  are  meant  to  reproduce 
well  only  the  ground  state  energy  and  the  first  excited  state  of  the  alkali  atom. 

As  a first  numerical  test  of  the  accuracy  of  our  pseudopotential  we  present  some 
results  for  the  ionization  energies  of  Li  and  Na.  Table  2-3  shows  the  dependence 
of  the  ionization  energies  of  Li  on  the  basis  sets.  Basis  set  (4s4p/2s2p)  reproduces 
well  the  ground  state  energy  of  Li.  For  the  other  calculations  we  started  with  a 
contracted  (9s9p5d/7s7p5d)  basis  set3  and  the  remaining  choices  were  obtained  by 
removing  functions  from  that  original  basis. 

In  Table  2-4  we  summarize  our  best  ionization  energies  for  Li  and  Na.  Our 
results  are  in  excellent  agreement  with  the  experiment. 

2.3.4  Choice  of  the  Core-Core  Repulsion  Function 

Finally,  we  consider  the  alkali  ion-Rg  interaction,  which  is  necessary  to  obtain 
the  true  adiabatic  potentials.  The  short-range  term  V"(R)  occurring  in  Eq.  2.11 
is  the  most  difficult  one  for  evaluation  and  seems  to  be  the  main  source  of  error  in 
the  estimation  of  Ej(R ) in  the  intermediate  range  of  the  internuclear  separation.  It 
can  be  experimentally  deduced  from  mobility  and  beam  scattering  measurements. 
Theoretically  it  is  calculated  mainly  on  the  basis  of  the  electron-gas  model.  In  pre- 
vious calculations  it  has  been  obtained  on  the  basis  of  the  statistical  Gombas-Baylis 
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Table  2-3:  Li  atom  ionization  energies  (au)  for  different  basis  sets 


Term  Exp.59  7s6p3d  6s6p5d  5s5p4d  4s4p3d  2s2p 


Li(2s)  2S 
Li(2p)  2P 
Li(3s)  2 S' 
Li(3p)  2P 
Li(3d)  2D 
Li(4s)  2S 
Li(4p)  2P 


0.198142 

0.130236 

0.074182 

0.057236 

0.055606 

0.038616 

0.031975 


0.198030 

0.130197 

0.074086 

0.057169 

0.055548 

0.038514 

0.031895 


0.198029 

0.130197 

0.074086 

0.057168 

0.055548 

0.038510 

0.031891 


0.198026 

0.130197 

0.074053 

0.057162 

0.055493 

0.038517 

0.031652 


0.198023 

0.130197 

0.032669 

0.056870 

0.043098 

0.015093 

-0.14235 


0.198067 

0.130078 

0.040740 

-.0687640 


Table  2-4:  Alkali-atom  ionization  energies  (au) 


Li  term 

9s9p5d 

Exp.59 

Na  term 

10s8p4d 

Exp.59 

Li(2s)  2S 

0.198115 

0.198142 

Na(3s)  2S 

0.188837 

0.188858 

Li(2p)  2P 

0.130127 

0.130236 

Na(3p)  2P 

0.111559 

0.111560 

Li(3s)  2S 

0.074176 

0.074182 

Na(4s)  2S 

0.071579 

0.071579 

Li(3p)  2P 

0.057111 

0.057236 

Na(3d)  2D 

0.055786 

0.055941 

Li(3d)  2D 

0.055482 

0.055606 

Na(4p)  2P 

0.050944 

0.050938 

Li(4s)  2S 

0.038614 

0.038616 

Na(5s)  2S 

0.037483 

0.037585 

Li(4p)  2P 

0.031886 

0.031975 

Na(5p)  2P 

0.028008 

0.029199 
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model.34  In  more  recent  calculations  Kim  and  Gordon56  and  Patil57  reproduced  the 
minima  of  the  2 II  term  very  well. 

We  have  opted  instead  to  use  a pragmatic  approach  to  calculate  the  VCSJ(R) 
term:  first,  we  subtract  our  ground  state  PEC  of  the  alkali-Rg  system  (omitting 
the  energy  contribution  from  the  Vcc ) from  the  most  accurate  curve  available  in  the 
literature.  The  curve  we  obtain,  which  represents  the  core-core  contribution,  is  then 
fitted  using  a function  of  adjustable  parameters  A and  b 

Ap-bR 

Vcscr(R)  = (2.13) 

Here  R is  a scalar  that  stands  for  the  internuclear  distance  between  the  two  atoms. 
Figure  2-1  shows  how  the  parameters  were  obtained  for  the  Li-He  pair.  The  energy 
curve  2scr  without  Vcc  corresponds  to  the  molecular  ground  state  curve  of  the  Li-He 
without  the  core-core  contribution;  2scr  Behmenburg  is  the  ground  state  curve  for 
the  Li-He  system  taken  from  an  ab  initio  calculation  by  Behmenburg  et  al.;4  2 so 
with  Vcc  is  the  energy  curve  of  the  molecular  ground  state  including  the  core-core 
energy  contribution. 

2.4  Computational  Aspects 

2.4.1  Pseudopotential  Matrix  Elements  in  the  Static  Atomic  Basis 

The  Hamiltonian  matrix  elements  including  pseudopotential  (neglecting  spin- 
orbit  coupling)  are  obtained  projecting  the  Hamiltonian  operator  on  the  left  and 
right  by  (xP  \ and  | Xq)  respectively: 

n„  = (xp\hSp\x,) 

= (xP  I - + Va(ta)  + V„b(im,R)  I X,)- 


(2.14) 
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The  first  term  in  the  above  equation  is  the  kinetic  energy  whereas  the  second  term 
describes  the  electron-core  interaction  which  is  composed  of 

K = (Xp  I VA(rA)  \ Xq)  = (Xp  I ;r  + VjTV)  + V?(rA)  | *,).  (2.15) 

'A 

The  first  element  in  Eq.  2.15  gives  the  electron-core  coulombic  attraction.  The  third 
term  in  Eq.  2.14  is 

KB  = (Xp  I Vab(ta,  R)  | x,)  = (Xp  I Vb{tb)  + Vec(rA,  R)  + VCC(R)  \ Xq)  (2.16) 

The  short-range  matrix  elements  found  in  Eq.  2.15  and  also  in  the  Vb(tb ) term  of 
Eq.  2.16  (as  shown  in  Eq.  2.9)  can  be  written  as: 

( XP  I V£(rx)  | xq)  = (XP  I ^R/,iexP (-/Vx)pJ,x  I Xg),  (2.17) 

with  X = A,  B. 

2.4.2  Programming  Details 

The  molecular-structure  calculations  reported  in  the  next  section  have  been 
carried  out  using  a series  of  subroutines  we  coded  in  Fortran  90. 

Following  a recursive  algorithm  proposed  by  Obara  and  Saika60  we  coded  the 
routines  gkineticg. f 90  and  gnuclearg.f90  to  calculate  the  kinetic  energy  and 
electron-core  attraction  matrices  respectively. 

The  routine  short-range . f 90  based  on  a procedure  described  by  McMurchie 
and  Davidson61  with  some  additions  from  Pitzer  and  Winter62  allowed  us  to  calculate 
the  short-range  matrix  elements  in  Eqs.  2.15-2.16. 

The  polarization  terms  of  Eq.  2.15,  the  Vb(tb)  term  of  Eq.  2.16  and  the  cross 
terms  of  Vec(rA,  R)  were  calculated  with  the  polar,  f 90  routine.  Here  we  followed 
the  formalism  presented  by  Schwerdtfeger.63 

The  procedure  we  have  developed  to  calculate  the  adiabatic  potentials  is  part 
of  the  eiktdmo.f90  program.  See  appendix  A for  the  flow  diagram  of  the  program 
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and  for  a detailed  explanation  of  the  subroutines.  Complete  derivations  of  the 
integrals  in  the  polar. f 90  and  the  short_range.f90  are  found  in  appendices  B 
and  C respectively. 

Similar  routines  for  the  calculation  of  these  matrix  elements  are  used  by  other 
quantum  chemistry  packages  like  MOLPRO64  and  COLUMBUS.65 

2.5  Results  and  Discussion 

The  present  calculations  of  the  interatomic  potentials  were  performed  for  the 
internuclear  separation  interval  between  3.0  and  20.0  au  with  different  step  sizes. 
The  PECs  have  been  calculated  for  the  LiHe,  LiNe,  NaHe  and  NaNe  diatomic  sys- 
tems. To  test  the  validity  of  the  pseudopotential  method,  our  PECs  are  compared 
with  other  theoretical  and  experimental  results.  The  simplest  system  of  that  series, 
LiHe,  is  studied  in  detail  to  determine  the  effect  of  the  basis  sets  size  of  the  PEC. 
The  energy  contribution  of  the  polarization  terms  is  also  analyzed. 

2.5.1  Comparison  of  the  LiHe  Potential  Energy  Curves  for  Different 
Basis  Sets 

An  optimal  basis  set  is  one  that  minimizes  the  computational  effort  without 
sacrificing  the  quality  of  the  PECs  of  the  excited  molecular  states  of  interest.  To 
choose  a good  basis  set,  several  preliminary  calculations  with  different  basis  sets 
need  to  be  run  until  the  balance  effort-quality  is  reached.  As  an  example  we  have 
done  this  in  the  case  of  the  LiHe  system.  For  a detailed  comparison  we  have  chosen 
the  basis  sets: 

Basis  II  : (6s5p3d/4s4p3d). 

Basis  III  : (7s6p4d/5s5p4d). 

Basis  IV  : (9s9p5d/7s7p5d). 

Figures  2-2-2-4  show  the  PECs  of  LiHe  (for  the  states  up  to  3d  term)  obtained 
using  the  basis  sets  II,  III  and  IV  respectively.  The  general  shape  of  the  potential 
curves  can  be  characterized  as  follows.  All  states  except  2scr  and  the  2pcr  exhibit 
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pronounced  minima  in  the  vicinity  of  4 au.  The  2sa  and  the  2pcr  are  purely  repulsive. 
A close  comparison  of  figures  2-2-2-4  shows  that  basis  set  III  gives  converged  results 
up  to  the  3d-curves. 

In  figure  2-5  our  LiHe  PECs  (basis  set  IV)  and  those  from  Czuchaj  et  al.3  are 
compared.  Their  PECs  were  obtained  with  a /-dependent  pseudopotential  for  the  Li 
core,  while  the  He  electrons  were  described  explicitly.  The  agreement  of  the  curves 
is  good  except  in  the  region  between  3 and  4 au  where  the  repulsive  term  of  the 
Hamiltonian  plays  an  important  role. 

Figure  2-6  compares  our  results  with  basis  set  IV  with  the  most  current  ab 
initio  results  from  Behmenburg  et  al.4  The  agreement  is  good  although  there  are 
some  deviations.  For  the  2p7 r state  our  De  = 1200  cm-1  and  Re  = 3.3  au  are  different 
from  the  ab  initio  values  De  = 868  cm-1  and  Re  = 3.42  au.  The  calculated  De  is 
compared  with  the  experimental  800  cm-1  by  Lee  et  al.66 

Later  in  the  chapter,  the  PECs  for  the  other  alkali-Rg  systems  will  be  pre- 
sented. 

2.5.2  Comparison  of  LiHe  Potential  Energy  Curves  by  Addition  of  Po- 
larization Terms  in  the  Hamiltonian 

We  now  turn  our  attention  to  the  polarization  terms  of  the  Hamiltonian  and 
their  impact  on  the  PECs.  All  calculations  were  performed  with  the  basis  set  IV  on 
Li  and  as  shown  in  figure  2-7,  where  all  the  PECs  have  been  labeled  as  either  (a), 
(b)  of  (c).  The  (a)  curves  were  obtained  with  the  full  Hamiltonian,  the  (b)  curves 
were  calculated  neglecting  the  polarization  cross-terms  and  the  (c)  were  produced 
neglecting  all  polarization  terms  and  cross-terms.  From  the  figure  we  conclude  that 
the  contribution  of  the  polarization  terms  to  the  PECs  is  appreciable  and  that  of 
the  cross-terms  negligible.  Therefore,  it  is  possible  to  neglect  the  cross-terms  from 
our  calculations  without  affecting  the  accuracy  of  our  results. 
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2.5.3  Calculation  of  the  Potential  Energy  Curves  of  the  LiNe,  NaHe  and 
NaNe  Systems  and  Comparisons  with  Other  Theoretical  Results 

Molecular  structure  calculations  have  been  performed  to  obtain  the  adiabatic 
potentials  for  the  ground  state  and  several  excited  states  of  the  LiNe,  NaHe  and 
NaNe  systems.  Internuclear  separations  ranging  from  3.0  to  25.0  au  with  different 
step  sizes  were  considered.  The  atomic  basis  sets  were  kept  constant  for  all  values  of 
distance  R.  The  precision  of  our  calculations  can  only  be  estimated  by  comparison 
with  spectroscopic  data  and  ab  initio  results. 

LiNe:  We  have  calculated  the  PECs  with  the  basis  set  IV  for  Li.  In  figure  2-8 
we  show  our  results  for  the  lower  states  of  Li  up  to  the  Li (3s)  term  along  with  those 
from  an  ab  initio  calculations  by  Behmenburg  et  al.31  The  agreement  of  the  curves 
is  good  since  our  potential  curves  present  minima  in  the  vicinity  of  4 au.  There 
are  some  minor  discrepancies  specially  in  our  2p<r  potential  curve  which  is  more 
repulsive  than  the  ab  initio  one. 

NaHe:  We  have  obtained  the  PECs  with  the  basis  set  (10sl0p3d/7s6p3d)  of 
Na.  Our  results  up  the  Na(4s)  state  are  compared  to  other  ab  initio  calculations  by 
Theodorakopoulos  and  Petsalakis67  in  Fig.  2-9.  Our  calculated  PECs  are  in  good 
agreement  with  the  ab  initio  ones  over  the  entire  distance  range,  except  our  3 pa 
curve  which  is  more  repulsive  . 

NaNe:  The  PECs  were  produced  with  the  same  basis  set  with  used  for  NaHe. 
In  figure  2-10  we  present  our  results  of  the  lowest  states  of  the  system  up  to  the 
3pcr  state.  Pseudopotential  results  by  Czuchaj  et  al.41  obtained  with  a basis  set 
(13sl2p4d/10s8p4d)  are  also  reported  in  that  figure.  There  were  no  ab  initio  results 
available.  The  agreement  of  the  curves  is  good  over  the  entire  range  of  distances. 
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Figure  2-1:  Vcc  curve  for  LiHe:  ^4=16. 016  and  b=  2.4402. 
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Figure  2-2:  LiHe  potential  energy  curves  using  basis  set  II. 
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Figure  2-3:  LiHe  potential  energy  curves  using  basis  set  III. 
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Figure  2-4:  LiHe  potential  energy  curves  using  basis  set  IV.  Note  that  it  is 
almost  identical  to  figure  2-3  within  the  graphics  accuracy. 
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Figure  2-5:  LiHe  potential  energy  curves  comparison:  Our  results  (basis  set 
IV)  vs.  pseudopotential  results  from  Czuchaj. 
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Figure  2-6:  LiHe  potential  energy  curves  comparison:  Our  results  (basis  set 
IV)  vs.  ab  initio  from  Behmenburg. 
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Figure  2-7:  Comparison  of  the  potential  energy  curves  of  LiHe.  (a)  results 
obtained  keeping  all  the  terms  of  the  Hamiltonian;  (b)  results  ob- 
tained neglecting  the  cross-terms;  (c)  results  obtained  neglecting 
all  the  polarization  terms. 
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Figure  2-8:  LiNe  potential  energy  curves  comparison:  Our  results  (basis  set 
IV)  vs.  ab  initio  from  Behmenburg. 
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Figure  2-9:  NaHe  potential  energy  curves  comparison:  Our  results  vs.  ab 
initio  results  from  Theodorakopoulos. 
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Figure  2-10:  NaNe  potential  energy  curves  comparison:  Our  results  vs. 
pseudopotential  results  from  Czuchaj. 


CHAPTER  3 

DYNAMICS  OF  ALKALI  ATOM-RARE  GAS  ATOM  COLLISIONS  IN  THE 

KEV  RANGE 

3.1  Introduction 

We  introduce  here  a new  first  principles  treatment  to  study  the  time  evolution 
of  electronic  transitions  in  collisions  of  alkali  atoms  with  Rg  atoms.  Our  approach 
combines  the  /-dependent  pseudopotential  method  (described  in  the  previous  chap- 
ter) and  the  eikonal/time  dependent  molecular  orbital  (Eik/TDMO)  method24  to 
reduce  our  many  electron  alkali-Rg  atom  problem  into  a three-body  problem  (alkali- 
metal  core,  valence  electron,  and  Rg  atom),  simplifying  dramatically  the  theory  and 
decreasing  the  computational  effort. 

We  begin  this  chapter  by  showing  the  derivation  of  the  Eik/TDMO  method  in 
terms  of  /-dependent  pseudopotentials  (Eik/TDMO-PP).  We  then  continue  present- 
ing a general  description  of  the  computational  aspects  of  the  method  and  a brief 
explanation  of  the  program.  It  is  of  interest  to  describe  what  is  occurring  during 
the  collisions.  Among  the  possible  quantities  of  interest  we  present  in  this  chapter 
results  for  the  orbital  populations  and  the  alignment  and  orientation  parameters 
of  the  collisional  system.  Later,  to  check  the  validity  of  our  method,  we  compare 
our  calculated  integral  cross  sections  (ICSs)  with  other  experimental  and  theoretical 
results. 

3.2  Eikonal/Time  Dependent  Molecular  Orbital  Method  Including 

/-Dependent  Pseudopotentials 

The  Eik/TDMO  method  is  the  version  of  the  more  general  eikonal/time  de- 
pendent Hartree-Fock  (Eik/TDHF)  method  to  study  one-active-electron  systems. 
The  implementation  of  Eik/TDHF  and  Eik/TDMO  methods  has  been  described  in 
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great  detail  in  references23,24,68  and  therefore,  we  only  summarize  the  most  relevant 
equations  of  the  Eik/TDMO  method  and  present  them  in  the  context  of  /-dependent 
pseudopotentials. 

We  consider  the  collision  of  an  alkali  atom  A with  an  Rg  atom  B.  The  total 
wave  function  for  the  system  is  a solution  of  the  Schrodinger  equation  for  total 
energy  E 


H*(X,R)  = EV(X,R), 


(3.1) 


where  X is  the  collection  of  coordinates  of  the  active  electron  and  R the  relative 
position  of  the  atomic  centers. 

The  Hamiltonian  operator  H of  the  atomic  pair  can  be  written  as 


H — — 


h2  d2 
2 M dR? 


(3.2) 


where  M is  the  reduced  mass  of  the  nuclei  and  H ^ is  the  Hamiltonian  term  for  fixed 
nuclear  positions.  In  what  follows,  we  have  replaced  by  the  Hamiltonian  Hpp 
(Hamiltonian  including  the  pseudopotential). 

It  is  convenient  to  write  the  total  wave  function  in  the  Eikonal  representation 


*(X,i2)  = x(X,J2)  exp  [iS{R)/h], 


(3.3) 


where  S is  real  mechanical  action  and  x is  the  complex  electronic  wave  function. 
Replacing  ^ in  Eq.  3.1,  we  obtain 


ypp 


E 


X(R))  = 0. 


2 M\idR  dR , 

Taking  the  real  part  of  the  projection  of  x on  the  right  of  Eq.  3.4  gives 


(3.4) 


2M  \dR 


OS 


(3.5) 
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where  Vqu  is  a quantum  potential 


Vqu  = v + v'  + v", 

V = 


<x  I A"  I x) 


(x  I x) 


V'  = 


V"  = 


ih  dS  / d\ 
2M  dR  |.\<9i? 

/ d2x 


x - u 


h2  1 


dx 

dR/ J 

^2x\ 

dR2/. 


(x  I x)> 
(x  I x)- 


(3.6) 

(3.7) 

(3.8) 

(3.9) 


2M2  [\dR 2 

The  first  term  of  the  quantum  potential  V^u  is  identified  as  the  Ehrenfest  po- 
tential for  core  motions,  while  the  second  and  third  terms  give  quantal  corrections 
due  to  nuclear  displacements. 

The  trajectories  of  the  nuclei  can  be  obtained  for  a fixed  energy  E introducing 

position  and  momentum  functions  of  the  time  t,  satisfying  the  Hamiltonian  equations 

— * —* 

P-P 


Hqu(P,R ) 


2 M 

dR  dH, 
dt 
dP_ 
dt 


+ Vqu{P , R )> 


qu 


dP 


dH, 


qu 


dR 


(3.10) 

(3.11) 

(3.12) 


where  we  defined  the  momentum 


P = 


dS_ 
dR ’ 


(3.13) 


and  the  total  mechanical  action  along  the  trajectory 

S(R)  = S(Ri)+  [ dR  ■ P. 

JRi 


(3.14) 


— r 

The  eikonal  approximation  for  short  de  Broglie  wavelengths  A = h/  \ Pr  \ is 
introduced  by  neglecting  the  second  and  third  terms  in  the  quantum  potential  that 
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involve  the  gradient  and  Laplacian  of  the  electronic  wave  functions: 

|)  / I (X  I X)  l«  A-*, 


(3.15) 


and 


(x 

d2x\ 

/ 

L 

dx  \ 

V 

OR2/ 

/ 

V 

OR2/ 

« A 


-l 


(316) 


which  allows  us  to  neglect  V'  and  V"  so  that  Vqu  « V . Hence,  we  propagate  the 
nuclei  under  the  influence  of  the  Ehrenfest  potential. 

Returning  to  the  description  of  the  electronic  states,  the  time-dependent  elec- 
tronic wave  function  for  our  one-electron  system  x is  reduced  to  a single  (complex 
valued)  molecular  orbital  ifj.  The  time  evolution  of  the  molecular  orbitals  can  be 
obtained  by  either  solving  a differential  equation  in  time  in  terms  of  the  molecular 
orbitals  ip  or  equivalently  by  solving  one  in  terms  of  the  electronic  density  operator 


p(t)  =i  ip(t))  i . 


(317) 


The  collision  dynamics  is  carried  out  solving  the  Hamiltonian  equations  3.11- 
3.12  coupled  to  time-dependent  differential  equation  for  the  density  operator  p 

- pH™  = ihBp/at.  (3.18) 


The  quantum  potential  in  terms  of  the  density  matrix  is  therefore  approximated 


by 


tr  [fiHtn 
tr[p] 


(3.19) 


and  becomes  a function  of  time. 


37 


3.2.1  Matrix  Equations 

The  time  evolution  of  electronic  densities  and  nuclear  positions  and  momenta, 
within  the  framework  of  the  Eik/TDMO-PP  method,  requires  the  choice  of  a basis 
set.  We  expand  the  time-dependent  molecular  orbitals  ip  in  terms  of  linear  combi- 
nations of  static  atomic  functions  (AFs)  {x},  which  are  centered  on  the  alkali  atom. 
It  is  well  known  that,  by  replacing  these  AFs  in  Eq.  3.18,  spurious  asymptotic  cou- 
plings appear  in  the  form  of  ( x p | dxp'/dt)  terms.23  These  spurious  couplings  are 
avoided  by  expanding  each  TDMO  ip  as  a linear  combination  of  the  traveling  atomic 
functions  £: 

t)  = ^2Ur,t)cip(t),  (3.20) 

p 

M?,t)  = Xp(r)Tm(r,  t),  (3.21) 


here  the  c’s  are  complex  expansion  coefficients,  Xp(r)  is  an  atomic  orbital  centered 
at  core  position  Rm(t ) for  the  electron  with  position  r.  The  electron  translation 
factor  (ETF)  Tm(r,t ) is  defined  by 


Tm(r,t)  = exp 


imP 


.(*) 


J'dtvKt')/ 2) 


(3.22) 


where  vm  is  the  velocity  vector  of  core  m,  ti  is  the  initial  time,  and  me  is  the  electron 
mass.  In  what  follows  we  are  setting  h = 1. 

The  density  operator  in  the  basis  of  traveling  atomic  orbitals  is  written  as 


m = Y,  I (p)PM(U  1=1  «>P«  I,  (3.23) 

pg 

where  the  P is  the  density  matrix.  The  matrix  elements  of  P are 


Rpq(t)  — ^ 'Cpj{t)Cqi(t)- 


OCC  l 


(3.24) 
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Equation  3.18  can  be  written  in  terms  of  the  density  matrix  by  projecting  on 
the  left  and  right  by  (£  | and  | £)  respectively  and  gives 

iP  = S-xHtP  - PHrS-1.  (3.25) 

Here  S_1  is  the  inverse  of  the  overlap  matrix  and  in  more  detail,  (HT)^  is  written 
as 

(H T)„  = ((„  | TmHSFT~'  | („),  (3.26) 

where  Tm  is  the  ETF  on  center  m.  It  is  now  clear  that  we  have  to  calculate  the 
one-electron  integrals,  present  in  the  matrix  form  of  Eq.  3.18,  in  a basis  of  traveling 
atomic  functions. 

We  discuss  now  the  transformation  of  matrices  from  a static  basis  \ 0n  which 
they  are  originally  calculated)  to  the  traveling  basis  £.  We  define  the  transformation 
matrix 


' = E (S(x,)m(S<x«)».. 

b 

(3.27) 

K 

(Sw>)„,,  = {X,  1 (,), 

(3.28) 

so  that  the  relationship  between  a matrix  in  the  static  basis  and  one  in  the 
traveling  basis  is  given  by 

0($)  = (3.29) 


3.2.2  Propagation  of  the  Density  Matrix 

The  propagation  of  the  density  matrix  is  carried  out  using  the  relax-and-drive 
procedure.23,68  Let  us  begin  by  rewriting  equation  3.25  in  the  form 


iP(t)  = W(t)P(<)  - [P(()W(()]t, 


(3.30) 
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where  W = S-1(f)H(t). 

The  density  matrix  P is  then  separated  into  two  matrices:  The  first  matrix, 
P(°\  describes  the  electronic  relaxation  for  fixed  nuclei,  and  the  second  matrix,  Q, 
describes  the  density-matrix  change  due  to  the  driving  forces  of  the  nuclear  motions: 

P(t)  = P<°>(t)  + Q(t).  (3.31) 

The  reference  density  P(°)  is  defined  by  its  time  propagation 

iP(0)(0  = W(t0)P(0)(t)  - P(0)(t)[W(t0)]t,  (3.32) 

or  equivalently, 

P°(f)  = U0(t,  t0)P0(f0)Uo(t,  to)1,  (3.33) 

where  U0  is  the  time  evolution  operator  defined  as 

U0(t,  t0)  = exp[— i(t  - f0)Wo].  (3.34) 

The  reference  density  at  the  beginning  of  a time  interval  is  given  by 

P(to)  = P(0)(to)  + Q(to).  (3.35) 

Given  this  choice  of  the  reference  density,  the  time  propagation  of  the  density-matrix 
change  Q follows  from 

*Q(<)  = W(t0)Q(t)  - Q(t)[W(t0)]t  + D(t),  (3.36) 

where  D is  the  driving  term 


D(t)  = AW(t)P°(t)  - P°(t)[AW(t)]t, 
AW(t)=W(t)-W0(to), 


(3.37) 

(3.38) 
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and  gives 


(3.39) 


The  propagation  of  the  above  differential  equations  requires  the  choice  of  a time 
step  A t = ti  — t0  which  is  chosen  so  that  the  condition 


Slower  ^ ||  P ||  ~ ^ higher  (3.40) 

is  satisfied.  The  lower  tolerance  tiower  allows  for  larger  time  steps  in  the  adiabatic 
region,  where  the  evolution  of  the  density  matrix  is  nearly  equal  to  the  evolution 
of  the  reference  density.  The  higher  tolerance  e/,j9/,er,  restricts  the  procedure  to 
smaller  time  steps  in  the  diabatic  region  where  the  change  in  the  density  matrix  is 
relatively  large.  The  relax-and-drive  procedure  is  able  to  adjust  the  time  step  size 
automatically  in  order  to  maximize  the  performance  of  the  method  by  minimizing 
the  number  of  steps  taken  for  the  maximum  accuracy  in  the  integration  the  density 
matrix  equations. 

3.2.3  Physical  Properties  of  Interest 

It  is  of  interest  to  describe  in  some  way  what  is  occurring  during  collisional 
electronic  excitation  and  rearrangement.  Among  the  possible  quantities  of  interest 
are  the  alignment  and  orientation  parameters  and  the  electronic  population  of  the 
atomic  orbitals. 

The  Eik/TDMO-PP  method  allows  us  to  propagate  the  density  matrix  over 
time  and  from  the  density  matrix  it  is  possible  to  examine  the  expectation  values  of 
an  operator  O by  using  the  general  expression 

tr[PQ] 
tr[P]  ‘ 


(6) 


In  this  study,  as  shown  in  Fig.  3-1,  we  consider  the  plane  of  the  collision  to  be  on 
the  x-z  plane  with  the  ^-direction  being  the  incoming  direction  of  the  projectile  and 
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the  x-direction  being  the  direction  of  the  impact  parameter,  so  that  the  y-direction 
is  always  perpendicular  to  the  plane. 

Having  chosen  this  particular  geometry,  we  move  on  to  the  definition  of  the 
alignment  and  orientation  parameters.  In  general,  one  may  describe  the  orbital  po- 
larization of  a system  by  examining  the  anisotropy  of  the  orbital  angular  momentum 
of  the  electrons.  A state  of  given  angular  momentum  L is  said  to  be:69 

1)  Isotropic,  if  all  its  magnetic  substates  are  equally  populated. 

2)  Oriented,  if  the  magnetic  substates  have  differing  populations. 

3)  Aligned,  if  substates  with  magnetic  quantum  numbers  of  equal  magnitude  but 
opposite  signs  are  equally  populated. 

To  illustrate  these  definitions  let  us  consider  a system  with  angular  momentum 
L.  We  introduce  the  spherical  irreducible  tensor  operators  L gn\  where  we  consider 
n = 1,2  and  —n  < q < n 


(3-42) 


where  we  define  in  Cartesian  coordinates 


i iLy. 


(3.43) 
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Let  us  define  now  the  orientation  vector  O ^ and  the  alignment  tensor  at 
time  t, 


v'TJlTT) 


Af\L,t)  = 


V6 


(3.44) 

(3.45) 


L(L  + l)"”'"5 

where  the  quantities  in  angle  brackets  are  the  expectation  values  which  can  be 
expressed  in  terms  of  the  density  matrix  and  hence  may  be  calculated  as  functions 
of  time  in  our  method 

tr[PL(">] 


(if).  = 


fa- [P) 


(346) 


In  our  geometry  the  only  alignment  parameter  with  nonzero  value  is 

if  = -L(3  l\  - i2). 


70 


(3.47) 

(3-48) 


The  time  evolution  of  the  atomic  orbital  population  can  be  obtained  from  the 
density  matrix  introducing  the  Lowdin  populations 


»i(‘)  = £ £[S(f)I/2WP(t)WS(f)1/2]Ap.  (3.49) 

k\ 

for  atomic  orbital  a. 

The  alignment  parameter  may  also  be  written  in  terms  of  Lowdin  populations 
as 


42)m  = «2p,(()-2n£,,(«).  (3.50) 

We  are  also  interested  in  reporting  final  physical  quantities  that  can  be  com- 
pared with  other  theories  and  experiments  such  as  state-to-state  integral  cross  sec- 
tions (ICSs).  The  ICS  from  an  initial  atomic  orbital  p to  a final  orbital  q for  projectile 
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energy  E is  defined  by 

ap^Q  (E)  = 2x[dbbpW(E,b)  (3.51) 

where  Pq%\E,b)  is  the  transition  probability  for  a transition  p — > q and  b is  the 
impact  parameter. 

3.3  Computational  Aspects 

3.3.1  Programming  Details 

The  pseudopotential  code  was  entirely  written  in  Fortran  90  (F90).  To  make 
possible  the  further  integration  of  the  pseudopotential  routines  into  the  original 
eiktdmo  code,  available  in  our  group,  it  was  necessary  to  rewrite  many  of  the  original 
Fortran  77  routines  into  Fortran  90. 

We  present  in  appendix  A a flow  diagram  for  the  algorithm  of  the  eiktdmopp 
program  along  with  an  explanation  of  the  main  features  of  each  subroutine. 

3.3.2  Construction  of  the  Density  Matrix 

The  matrix  representation  of  the  density,  overlap  and  Hamiltonian  operators 
requires  a set  of  basis  functions  in  which  the  atomic  and  molecular  orbitals  for  the 
ground  and  several  other  excited  states  may  be  expanded  accurately. 

Our  TDMOs  {^}  are  written  as  combinations  of  atomic  Cartesian  functions 

m 

I i>)  = I <PC)CC,  (3.52) 

where  Cc  is  the  matrix  of  the  expansion  coefficients.  Here  we  have  assumed  that 
the  atomic  functions  are  normalized  but  not  orthogonal  to  each  other.  The  matrix 
Cc  is  the  matrix  of  eigenvectors  obtained  by  diagonalizing  the  Wc  matrix  (in  the 
{4>c}  basis)  at  initial  time 


WCCC 


CCE. 


(3.53) 
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Once  we  have  selected  the  initial  state  of  the  system  ipi  (with  i = 2s,  2 px, 
etc.),  which  in  our  calculations  corresponds  the  ground  state  of  the  alkali  atom,  we 
construct  the  initial  density  matrix  (in  the  basis  of  molecular  orbitals  ip)  setting 
all  its  matrix  elements  to  zero  except  the  element  which  is  set  to  one.  Because 
all  the  matrices  we  use  in  the  propagation  are  calculated  in  the  basis  of  Cartesian 
Gaussian  functions  (pc,  we  need  to  take  the  density  matrix  we  have  constructed  in 
the  ip  basis  to  the  <pc  basis  in  which  the  calculation  are  carried  out.  The  procedure 
we  have  developed  to  do  this  transformation  is  discussed  below. 

The  density  operator  p can  be  expanded  in  terms  of  either  the  TDMOs  ip  or 
the  basis  of  atomic  Cartesian  functions  (pc  (the  basis  in  which  the  calculations  are 
run) 

P = | | = | ^C)P‘(#  | . (3.54) 

We  transform  the  density  matrix  from  the  ip  basis  to  <pc  basis  by  projecting  the 
second  and  third  terms  of  Eq.  3.54  on  the  left  by  ( <pc  | and  on  the  right  by  | (pc ) to 
get 

(P  I <A'>Pc(^  I tf>c)  = (<f  I i>)P*(4>  I <t>c)  (3.55) 

SCPCSC  = (4f  | tWP*{®  | tf)  (3.56) 

and  replacing  the  ip’s  in  3.56  using  their  expansions  shown  in  3.52  we  get 

gcpcgc  = j ^Ccp*Cct^c  | ^ 

= scCcPv'CctSc.  (3.57) 

Multiplying  on  the  left  and  right  of  Eq.  3.57  by  (Sc)-1  we  get  the  expression  we 
were  looking  for  Pc  in  terms  of  P^ 


pc  _ 0cpV'(2;ct 


(3.58) 
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where  Pc  is  the  density  matrix  in  the  Cartesian  Gaussian  basis. 

Once  the  propagation  has  ended,  we  take  the  density  matrix  from  the  (pc  basis 
to  the  ip  basis  to  analyze  the  results,  using  the  transformation 

P*  = (Cc)-1Pc(Cct)-1.  (3.59) 

3.4  Results  and  Discussion 

In  this  section  we  report  some  of  the  results  we  have  obtained  by  applying 
our  method  to  the  LiHe,  LiNe,  NaHe  and  NaNe  collisional  systems.  Energies  of 
interest  for  the  incident  alkali  atom  range  from  as  low  as  0.5  keV  to  greater  than 
30  keV  depending  upon  the  system.  The  simplest  of  the  alkali-Rg  system  consisting 
of  Li  and  He  has  been  used  to  study  in  detail  the  time-dependence  of  the  orbital 
population,  orientation  and  alignment  parameters.  Later,  we  summarize  the  results 
of  final  excitation  ICSs. 

A careful  selection  of  the  initial  conditions  guarantees  the  convergence  and 
reproducibility  of  the  final  results.  Series  of  test  runs  ensures  the  convergence  and 
stability  of  the  numerical  procedures  for  the  integration  of  the  equations.  These 
tests  involve  varying  the  initial  and  final  separation  of  the  cores,  the  lower  and 
higher  tolerances.  Later,  in  section  3.4.2  we  study  the  effects  of  the  basis  set  size  on 
the  ICSs,  where  the  basis  set  choice  plays  an  important  role  in  the  convergence  of 
the  results. 

3.4.1  Time-Dependent  Study  of  the  Li-He  Collision 

As  discussed  in  chapter  2,  the  computational  effort  taken  in  the  calculation 
of  Hamiltonian  matrix  rises  rapidly  with  the  number  of  atomic  orbitals  described 
by  the  basis  set.  Therefore,  we  have  to  find  an  optimal  basis  set  that  reduces  the 
computational  time  without  sacrificing  the  reproducibility  and  convergence  of  the 
results.  There  are  other  initial  parameters  that  also  affect  the  final  results  of  the 
calculations  and  they  need  to  be  investigated  as  well.  We  have  found  that  initial  and 
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final  distances  of  30  au  (where  the  interaction  of  the  two  atoms  is  negligible)  are  good 
for  most  systems.  The  tolerance  set  plays  an  important  role  in  final  outcome.  It  has 
been  found  that  the  pair  10~7  — 10-5  (for  lower  and  higher  tolerances  respectively) 
gives  accurate  results  in  most  cases. 

We  now  report  some  of  the  results  we  have  obtained  from  studies  on  LiHe71  : 

Figure  3-2  presents  the  trajectories  of  the  nuclei  at  different  impact  parameters 
(6  = 0.5,  0.75,  1.0  au)  for  a fixed  initial  laboratory  energy  ( Eiab ) of  1.0  keV.  It  is 
important  to  note  in  Figure  3-2  how  the  He  atom  (initially  at  rest)  recoils.  We  also 
find  that  the  deflection  angles  are  small  even  for  these  small  impact  parameters.  This 
is  due  to  the  high  collisional  energy  we  have  chosen  and  also  the  neutral  character 
of  the  colliding  atoms  that  makes  the  core-core  repulsion  significant  only  when  the 
two  nuclei  overlap  at  short  distances. 

Figure  3-3  shows  the  population  of  the  Li(2s)  and  Li(2p)  orbitals  as  a function 
of  time  for  and  Eiab  of  1 keV.  It  is  always  instructive  to  analyze  the  temporal  change 
of  the  orbital  populations,  to  learn  about  the  nature  of  electronic  rearrangement. 
All  populations  oscillate  over  time  before  coming  to  a constant  asymptotic  value 
that  depends  on  the  energy  of  collision  and  the  impact  parameter.  The  population 
of  the  2 py  orbitals  is  not  shown  as  it  is  never  populated  by  the  chosen  geometry  of 
the  calculation.  The  total  time  intervals  over  which  populations  change  are  about 
200.0  au. 

Figure  3-4  gives  the  dependence  of  the  final  orbital  population  on  the  impact 
parameter  at  a Eiab  of  1 keV.  We  see  how  the  excitation  from  the  2s  to  the  2p 
orbitals  is  significant  in  the  region  of  impact  parameters  from  0. 0-6.0  au.  We  have 
found  that  this  feature  is  repeated  at  other  collisional  energies. 

Other  time-dependent  phenomena  that  are  of  interest  for  this  system  include 
the  alignment  parameter  presented  in  figures  3-5-3-7  and  orientation  parameter 
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presented  in  figures  3-8-3-10.  These  results  show  that  both  the  alignment  and  the 
orientation  parameters  also  oscillates  during  the  collision. 

3.4.2  Basis  Set  Effects 

The  validity  of  our  method  in  the  study  of  the  electronic  rearrangement  in 
alkali-Rg  collisions  is  tested  when  we  compare  our  ICSs  results  with  those  from 
other  theories  and  experiments. 

We  turn  our  attention  now  to  the  calculation  of  the  2s  — > 2 p integral  cross 
section  for  Li-He  system.  Calculations  were  performed  with  kinetic  energies  Eiab 
chosen  to  be  equal  to  experimental  values,  and  impact  parameters  in  the  range 
b — 0.0  — 20.0  au,  with  eighty  values  per  energy,  and  a higher  density  grid  at  lower 
impact  parameters.  Integrals  over  b converged  at  an  upper  limit  of  12.0  au.  For  the 
detailed  calculations  we  have  chosen  the  basis  sets: 

Basis  I:(6s5p2d/4s4p2d). 

Basis  II:(6s5p3d/4s4p3d). 

Basis  III:(7s6p4d/5s5p4d). 

Basis  IV:(9s9p5d/7s7p5d). 

Figure  3-11  gives  our  integral  cross  sections  for  Li-He  in  the  range  Eiab  = 1.0- 
10.0  keV,  comparing  basis  sets  I-IV  to  experimental  data.14  We  find  that  there  is  a 
strong  dependence  of  the  ICSs  on  the  size  of  the  basis  set.  The  sequences  of  larger 
basis  sets  are  found  to  converge.  This  comparison  shows  that  our  results  with  basis 
sets  III  and  IV  are  in  excellent  agreement  with  experiment.  The  experimental  values, 
which  are  not  measured  in  absolute  terms,  have  been  adjusted  to  our  calculations 
at  a single  energy  of  Eiab  = 5.0  keV.  These  tests  calculations  have  shown  that  to 
correctly  obtain  the  2s  — > 2 p ICS  it  is  necessary  to  include  several  d-orbitals. 
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Figure  3-12  shows  our  calculated  ICSs  (for  basis  set  III)  and  compares  them 
with  results  from  calculations  using  close-coupling  equations  for  a trajectory  of  con- 
stant velocity,  and  either  a model  (Baylis)  core  potential72  , or  atomic  pseudopo- 
tential.35 The  curves  show  the  same  general  behavior  with  energy  and  the  same 
orders  of  magnitude  of  ICSs,  but  noticeable  differences.  Our  results  also  allow  for 
close-coupling  and  use  a calculated  trajectory,  but  the  different  trajectories  are  not 
likely  to  change  ICSs  much,  so  that  the  likely  sources  of  differences  are  the  choice 
of  core  potentials  and  the  size  of  atomic  basis  sets. 

3.4.3  Integral  Cross  Sections  for  the  Li-Ne,  Na-He  and  Na-Ne  Systems 

We  now  move  on  to  ICSs  for  the  Li-Ne,  Na-He  and  Na-Ne  collisional  systems 
for  energies  in  the  keV  range  and  compare  them  with  experimental  data  and  results 
from  other  theories. 

Li-Ne:  Figure  3-13  shows  a comparison  of  our  ICSs  with  experimental  data14 
and  theoretical  results.72  The  experimental  values  have  been  adjusted  to  our  cal- 
culations at  the  single  energy  of  Eiab  = 3.0  keV.  Note  that  the  agreement  between 
calculated  and  experimental  results  is  good  through  the  entire  range  of  energies.  Our 
ICSs  are  also  compared  with  other  theoretical  ICSs.  We  find  that  the  curves  have 
the  same  general  behavior  with  energy.  However,  there  are  pronounced  differences 
in  magnitude.  These  differences  may  be  attributed  to  the  choice  of  model  potential 
in  those  calculations. 

Na-He:  This  system  has  been  the  subject  of  extensive  experimental  and  theo- 
retical investigations.  Figure  3-14  presents  the  ICSs  we  have  obtained  with  the  Na 
(10sl0p3d/7s6p3d)  basis  set.  We  include  experimental14  and  theoretical35  results. 
We  adjusted  the  experimental  results  to  our  calculation  at  the  single  energy  of  Eiab  = 
15.0  keV.  This  comparison  shows  that  our  theoretical  results  are  in  good  agreement 
with  both  experiment  and  theory  for  the  entire  energy  range. 
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Figure  3-1:  The  collision  is  taken  to  occur  in  the  x — z plane,  the  alkali  atom 
is  chosen  to  be  the  projectile  and  we  assume  that  the  rare-gas 
atom  is  initially  at  rest. 

Na-Ne:  This  is  the  last  of  the  systems  we  have  studied.  Figure  3-15  gives  the 
ICSs  obtained  with  the  Na  (10sl0p3d/7s6p3d)  basis  set.  Our  results  are  compared 
with  experimental  data  from  Olsen14  that  we  have  adjusted  at  the  energy  of  Eiab 
= 10.0  keV.  Both  curves  present  the  same  behavior  at  these  energies  although  it  is 
clear  that  the  trends  of  the  two  curves  at  energies  greater  than  15  keV  are  different. 
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Figure  3-2:  Trajectories  for  a 1.0  keV  Li  atom  incident  on  a He  atom  target. 

Comparison  of  results  for  three  different  impact  parameters,  b — 
0.5,  0.75  and  1.0  au.  The  He  target  is  initially  sitting  at  rest  at 
the  origin  of  coordinates. 
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Figure  3-3:  Populations  of  the  Li  2s  and  2p  orbitals  vs.  time  in  a 1.0  keV 
Li-He  collision,  at  b = 1.0  au. 
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Figure  3-4:  Final  orbital  populations  of  Li  2s  and  2p  orbitals  vs.  impact 
parameter  in  1.0  keV  Li-He  collisions. 
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Figure  3-5:  Alignment  parameter  vs.  time  in  a Li-He  collision  ( Eiab  — 1.0 
keV,  b = 0.5  au). 
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Figure  3-6:  Alignment  parameter  vs.  time  in  a Li-He  collision  ( E\ab  = 1.0 
keV,  b = 1.0  au). 
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Figure  3-7:  Alignment  parameter  vs.  time  in  a Li-He  collision  ( Eiab  = 1.0 
keV,  b = 2.0  au). 
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Figure  3-8:  Orientation  parameter  vs.  time  in  a Li-He  collision  ( Eiab  = 1.0 
keV,  b = 0.5  au). 
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Figure  3-9:  Orientation  parameter  vs.  time  in  a Li-He  collision  (Eiab  = 1.0 
keV,  b = 1.0  au). 
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Figure  3-10:  Orientation  parameter  vs.  time  in  a Li-He  collision  (Eiab  — 1.0 
keV,  b = 2.0  au). 
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Figure  3-11:  2s  — >2 p Integral  cross  sections  vs.  Eiab  for  Li-He  collisions. 

Comparison  of  our  results  using  different  basis  set  choices  with 
experiment. 
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Figure  3-12:  2s  —>■2 p Integral  cross  sections  for  collisions  of  Li-He  (Eiab  = 
1.0-10.0  keV).  Our  results  (using  basis  set  III)  compared  with 
other  theories. 


ICS  (A2) 
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Figure  3-13:  2s  — >2p  Integral  cross  sections  in  Li-Ne  collisions.  Compar- 
ison of  our  results  with  experimental  results  from  Olsen  and 
theoretical  results  from  Manique. 
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Figure  3-14:  3s  — >3 p Integral  cross  sections  in  collisions  of  Na-He.  We 
present  our  results  and  those  from  Olsen  (exp.)  and 
Kimura(theory). 


ics  (A2 ) 
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Figure  3-15:  3s  — >-3 p Integral  cross  sections  for  Na-Ne  collisions.  We  show 
our  theoretical  results  and  those  from  Olsen  (exp.). 


CHAPTER  4 

ALKALI  ATOMS  FINE-STRUCTURE  TRANSITIONS  INDUCED  BY 
COLLISIONS  WITH  RARE  GAS  ATOMS 

4.1  Introduction 

In  this  chapter  we  are  interested  in  studying  the  time  dependence  of  the  thermal 
energy  collisions  of  alkali  atom  with  rare-gas  atom  to  gain  insight  into  the  mechanism 
of  the  fine-structure  transitions  among  fine-structure  levels  of  the  alkali  atom. 

This  chapter  starts  with  a general  review  of  the  literature  since  this  is  the 
subject  of  numerous  experimental  and  theoretical  investigations.  Then,  the  theory 
of  spin-orbit  coupling  in  atoms  is  discussed  and  extended  to  diatomic  interactions. 
Some  computational  aspects  are  discussed  including  the  derivation  of  the  spin-orbit 
coupling  matrix  elements  in  the  static  basis.  Some  preliminary  results  of  our  calcula- 
tions including  spin-orbit  coupling  are  presented.  We  then  present  the  Eik/TDMO- 
PP  method  and  extend  it  to  include  spin-orbit  coupling  (Eik/TDMO-PP-SO).  A 
quick  overview  of  the  transformation  to  atomic  eigenstates  with  spin-orbit  coupling 
is  presented,  followed  by  a discussion  about  the  molecular  spin-orbitals.  Finally, 
some  aspects  related  to  application  to  the  Na*  -I-  He  system  are  discussed. 

4.2  Background 

The  thermal  and  hyperthermal  collisions  of  excited  alkali  atoms  M*  in  a sub- 
state 2Pj  with  Rg  atoms  result  both  in  angular  momentum  reorientation 

m*CPj,Mj ) + Rg^  m*CPj,m'j ) + Rg  (4.i) 

and  in  angular  momentum  recoupling 

M*(2Pj ) + Rg  M*(2Pj ,)  + Rg  (4.2) 
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with  energy  transfer  A Ejj>  . These  are  the  only  processes  which  take  place  with  ap- 
preciable probability  in  the  range  of  energies  of  a few  eV11  where  spin-orbit  coupling 
plays  and  important  role. 

These  fine-structure  transitions  have  been  extensively  studied,  but  more  can  be 
learned,  especially  with  regard  alignment  and  orientation  effects  and  the  mechanism 
of  the  transfer,  using  a time-dependent  description.  One  of  the  simplest  cases  is  the 
process 

Na(Zp  3P3/2)  + He  -»  Na(3p  2P1/2)  + He.  (4.3) 

A complete  and  rigorous  theory  of  these  processes  involves  a calculation  of  the 
full  quantum  wave  function  ^(R,  r)  describing  the  motion  of  the  nuclei  and  active 
electron.  This  close-coupling  theory  begins  with  an  expansion  of  the  wave  function  in 
molecular  electronic  states,  and  ends  with  numerical  solution  of  coupled  equations  for 
nuclear  wave  functions  and  numerical  summation  over  coupled  angular  momentum 
states.  Accurate  cross  section  have  been  reported  by  Reid27  , Pascale  and  Olson,28 
Lemoine  et  al., 29  Schatz  et  al.30  and  Leo  et  al.73  using  this  approach.  However 
the  formulation  of  the  theory  and  the  long  numerical  codes  are  very  complex  and 
physical  insight  can  be  lost. 

A more  intuitive,  though  less  rigorous,  semiclassical  model  has  successfully  been 
used  to  study  these  systems.  Nikitin36  used  a strong  coupling  approximation  to  ob- 
tain analytic  formulas  for  the  fine-structure  transition  cross  sections  in  the  Na-Ar 
system.  Masnow-Seeuws,74  with  Roueff38  numerically  solved  more  exact  equations 
for  the  Na-He  system  using  a semiclassical  impact  parameter  method;  this  calcu- 
lation was  later  refined  by  Masnou-Seeuws  and  McCarroll39  to  take  into  account 
trajectory  effects.  More  recently  Kovalenko  et  al.40  used  a semiclassical  formalism 
to  monitor  the  expectation  values  of  the  orbital  and  spin  angular  momentum  vectors. 
Inspired  by  the  latter  we  will  try  to  study  the  mentioned  collisional  phenomena. 
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4.3  Spin-Orbit  Coupling  in  Atoms 

Spin  is  of  course  a relativistic  effect,  but  a classical  picture  nonetheless  does  help 
in  the  understanding  of  spin-orbit  coupling.  The  interaction  couples  together  the 
spin  (s)  and  the  orbital  (1)  angular  momentum  vectors,  and  reduces  the  degeneracy 
of  spectroscopic  terms  with  a given  / and  s,  but  a different  value  of  j , the  total 
angular  momentum.  The  interaction  of  the  electron’s  magnetic  moment  /is  with  a 
field  B is  given  by  fis  ■ B and,  in  the  absence  of  an  external  field,  B is  proportional 
of  the  vector  product  of  the  nuclear  field  E and  velocity  of  the  electron  v 


B = 


E x v 


(4.4) 


where  c is  the  speed  of  light.  Now  /is  = —(eh/2mc)cr  and  v = p /me,  where  me  is 
the  electron  mass,  so  that  the  interaction  will  be 


(eh/2m2c2)(T  - (Exp) 


(4.5) 


and  for  a central  field 

„ ldV 

E = — 7T- 

r or 

with  V being  the  potential. 

Writing  s = her / 2 and  1 = r x p for  one  electron,  the  expression  for  the  inter- 
action becomes 

(e2/m2c2)^^l-s.  (4.7) 

We  have  considered  the  nucleus  to  be  at  rest,  and  so  this  must  be  divided  by  2 
(the  Thomas  precession  factor). 

Since  spin  is  a relativistic  effect,  it  is  more  correct  to  start  with  Dirac’s  relativis- 
tic equation  for  the  electron.  The  reduction  of  the  Dirac  equation  to  non-relativistic 
form  and  the  use  of  atomic  units  (e  = H = me  — 1)  gives  the  spin-orbit  interaction 
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energy  for  a one-electron  atom  as75 


Hso 


^ldVl 

2 r dr 


s, 


(4.8) 


where  V is  the  Coulomb  potential  due  to  the  nucleus  and  a = 1/c  is  the  fine- 
structure  constant.  This  may  be  rewritten  as  £1  • s and  this  provides  a definition  of 
the  spin-orbit  coupling  constant  £ for  a one-electron  atom.  In  such  case  V = —Z/r 
and  so  dV/dr  — Z/r 2 

The  generalization  of  this  result  to  many  electron  atoms  is  far  from  straight- 
forward but  following  Condon  and  Shortley76  it  is  normal  to  reinterpret  V as  the 
screened  potential  due  to  the  nucleus  and  to  all  the  other  electrons  and  the  result 
is  then  summed  over  all  the  outer  electrons,  giving  the  spin-orbit  Hamiltonian  as 


Hso  = 


2 i \ dri  1 Sl 


(4.9) 


This  spin-orbit  Hamiltonian  does  not  commute  with  the  operators  Sz,  Lz,  but 
with  S2,  L2,  J 2 and  Jz.  The  latter  pair  are  related  to  the  others  by 


Jz  — Lz  + Sz 


(4,10) 


and  its  analogues  and  by 

J2  = J2x  + J2y+J2z.  (4.11) 

If  we  have  a function  ip  which  is  an  eigenfunction  of  these  operators,  then 


Jzip  = Mjip 


Mj  — Mi  + Ms 

J2ip  = J(J  + l)iP 


(4.12) 

(4.13) 

(4.14) 

(4.15) 


-J  < Mj  < J . 
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In  order  to  calculate  matrix  elements  of  this  Hamiltonian  including  Hso  it  is  con- 
venient to  use  eigenfunctions  of  J2. 

4.4  Spin-Orbit  Coupling  in  Diatomic  Molecules 
There  is  no  rigorous  first-principle  derivation  of  the  spin-orbit  coupling  for  the 
two-center  many-electron  case.  Usually,  it  has  been  constructed  from  the  additivity 
of  the  Breit-Pauli  two-particle  Hamiltonian.75,77  For  a diatomic  molecule  one  has 
that 


where  c is  the  speed  of  light,  lj*  = r ix  xp;,  p*  is  the  momentum  of  the  ith  electron, 
and  riX  is  position  of  the  ith  electron  respective  to  nuclei  X (X  — A,  B).  The 
first  two  terms  in  Eq.  4.16  are  one-electron  operators  and  represent  the  spin-orbit 
coupling  of  each  electron  in  the  field  of  the  two  bare  nuclei  with  charges  Za  and 
ZB , respectively.  The  third  term,  the  spin-other-orbit  interaction,  arises  from  inter- 
electronic  interactions  and  partially  counterbalances  the  field  of  the  bare  nuclei.  Its 
sign  is  opposite  to  that  of  the  one-electron  term. 

The  use  of  Eq.  4.16  gives  good  agreement  with  experimental  results  for  the 
fine-structure  splittings,75  but  can  not  be  used  within  our  one-electron  formalism. 
Therefore,  we  will  opt  for  a different,  more  convenient  approach  in  terms  of  pseu- 
dopotentials. 

4.4.1  Spin-Orbit  Coupling  in  Terms  of  Pseudopotentials 

The  potentials  obtained  directly  from  relativistic  atomic  wave  functions  have 
the  form 

oo  1+1/2 

VRBF  = E E U,fp(r)Ph,  (4.17) 

1=0  j=|I— 1/2| 

where  we  use  the  notation  from  Pitzer  and  Winter.62  REP  denotes  relativistic 
effective  potential,  and  the  Pij  are  (spin-dependent)  projection  operators.  The  P^ 
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are  more  difficult  to  work  with  than  the  spin-independent  projection  operators  Pi, 
but  to  relate  the  two,  a spin  operator  must  also  appear,78  and  it  does  so  in  the  form 
1 ■ s: 


jjREP 


OO 


oo 


U?REP(r)Pi  + Zi(r)Pi  i s Pt  = UAREP  + Hso, 

t=o  1=1 


(4.18) 


where  JJAREP  is  an  averaged  (with  2j  + 1 weighting)  relativistic  effective  poten- 
tial and  £i(r)  depends  on  the  difference  of  Ui,t+i/2(r)  and  Ui,i-\/2(r).  These  two 
terms  are  readily  identified79  as  core  potentials  and  spin-orbit  (potential)  operators 
respectively. 

The  core  potential  is  reduced  to  the  potential  described  in  chapter  1.  The  spin- 
orbit  operator,  is  reduced  to  a finite  sum  by  making  the  approximation  that  the 
/ • s are  independent  of  l for  l > L,  where  L is  one  more  than  the  highest  quantum 
number  of  the  core  orbitals  being  replaced.  With  the  use  of  the  completeness  relation 
property  of  the  Pi,  and  the  truncation  of  terms  higher  than  ^(r),  it  yields. 


L—l 

hso  = 6.M i • s + [6W  - «i(r)l p'  l spi’  <4-19) 

Z=1 

but  as  shown  in  previous  calculations80,81  the  approximate  form 

L—l 

pS°  = 'El(l(r)-SL(r)lPll-sfii  (4.20) 

1=1 

produces  accurate  results.  Pitzer  and  Winter78  have  proposed  a spin-orbit  operator 
of  the  form 


Hso 


v 

*max 

E 


2Ay»°(,-A) 

21  + 1 /,A 


1 • s Pi,, 


(4.21) 


with  1 and  s denoting  the  orbital  angular  momentum  and  spin  operators,  A the 
atomic  cores  and  Pi,\  are  spin-independent  projection  operator.  The  difference 
of  the  radial  parts  of  the  two-component  relativistic  pseudopotentials  Vt \PiPi/2 
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Table  4-1:  Parameters  (in  au)  of  the  spin-orbit  operators  for  Li  and  Na. 


l 

nk 

ABlk 

Pik 

Li 

1 

2 

0.8078 

0.000035 

2 

2.5500 

-0.000035 

1 

7.2535 

0.000869 

Na 

1 

2 

1.6068 

0.000086 

2 

22.5231 

0.023508 

1 

76.2365 

0.032039 

and  is  written  (similarly  to  the  radial  part  of  the  spin-orbit  averaged  pseu- 

dopotentials) in  terms  of  the  Gaussian  functions: 

AF,so(riA)  = £ exp(-&r?A).  (4.22) 

k 

We  have  omitted  the  atomic  core  index  A index  from  our  notation  since  the 
operator  is  only  applied  to  the  alkali  atom. 

Elimination  of  the  core  electrons  leads  to  a molecular  valence  Hamiltonian  (in 
au) 


+ v;fp  + hso. 


(4.23) 


4.4.2  Choice  of  Parameters 

The  exponential  parameters  /?/&,  the  linear  parameters  AB^.  and  the  exponents 
nk  from  Eq.  4.22  constitute  adjustable  parameters  which  are  optimized  in  a least- 
squares  sense  to  reproduce  reference  spin-orbit  splittings  derived  form  all-electron 
calculations.  Table  4-1  lists  the  SOC  parameters  for  Li  and  Na.  These  numbers 
were  taken  from  Pacios  and  Christiansen.80 
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4.5  The  Eik/TDMO  Method  with  Spin-Orbit  Coupling 

We  start  the  TDMO  approximation  by  expressing  the  spin-orbit  coupling  states 
rj\  with  A = ( Isjrrij ) as 


t7a(x,  t)  = Yl  l)  (4-24) 

j 

where  the  at"  are  taken  from  the  spin-orbit  coupling  states  at  initial  time  tin,  and 
Xj  = (r^,  Cj)  are  a collection  of  position  and  spin  variables. 

For  a given  initial  electronic  state  I and  trajectory  R (t),  one  has  that 

ipi(xi,t)  = (4-25) 

(V>7  I $j)  = SijS'Yj',  (4.26) 


where  the  ip]  are  TDMOs  for  electron  spin  7 = a,  /3.  In  preparation  for  the  inclusion 
of  spin-orbit  coupling  in  the  TDMO  approximation,  it  is  convenient  to  define  the 
electronic  density  operator  for  spin  projections  7, 7 ' — a, /3  as 

p77'w  =1  flmwfw  1 (4-27) 


satisfying  the  TDMO  equations  with  spin-orbit  coupling  (TDMO-SO) 

dp* 


[^77"p7"7'  - p77"#7"7']  = ih 


dt 


(4.28) 


with  H1  7 an  operator  from  the  one-electron  Hamiltonian. 

4.6  Eik/TDMO-SO  in  Terms  of  a General  Basis  Set 
In  this  section  a general  basis  set  of  electronic  orbitals  is  introduced  in  which 
the  MOs  are  expanded  as 


V>,7  (r,  t)  = ^(pP{v,t)clp(t) 


(4.29) 


with  the  coefficients  being  complex  valued. 
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In  what  follows,  a general  matrix  operator  for  spins  77',  O77'  forms  a 2x2 
supermatrix  O,  according  to  the  notation 


0 = 


^ Qaa  qq/3^ 


Q0a  Q00 


(4.30) 


The  p77  element  of  the  density  superoperator  can  be  written  in  the  basis  of  Eq. 
4.29  as 


pw(t) = E 1 


(4.31) 


PQ 


where  P77'  is  the  (pg)-element  of  the  density  matrix  P77\  Eq.  4.31  in  matrix 


PQ 

notation  is 


P77'  =|  07>P77W  I (4-32) 

where  | (fry)  and  {(fry  | are  row  and  column  matrices  with  spin-orbital  elements  </>P7, 
respectively. 

The  matrix  elements  of  P77,  are 

^'m  = E <*(*)<*«■•  (4-33) 

occ  i 

The  one-electron  Hamiltonian  matrix  H77'  is  given  by 

H77'  = (07  | H | 07').  (4.34) 

Alternatively,  the  one-electron  Hamiltonian  matrix  can  be  written  as 

H1’y  = (K  + V")^,  + Hjg  (4.35) 

with  K being  the  electron  kinetic  energy  matrix,  V^p  is  the  pseudopotential  matrix 
described  in  chapter  2,  and  HJq  the  spin-orbit  coupling  matrix. 
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In  order  to  derive  the  TDMO-SO  equations  for  the  density  matrix,  the  matrices 
Q and  S are  defined,  respectively,  by 

= (<h\jt4n')  = n5w,  (4.36) 

and 


S17"  = <07  | 0V>  = S<S,y.  (4.37) 

Starting  from  the  four  differential  equations  for  p77  in  Eq.  4.28,  and  multiplying 
them  by  (<£7  | from  the  left  and  by  | 4>y)  from  the  right,  the  following  four  differ- 
ential equations  are  obtained 

iF  = S-1(H  - *n,) P - P(H  - *n4)tS"1,  (4.38) 


where  the  supermatrices  are  defined  as 


and 


(HQQ  Hap\ 

RPa  H/3/3  J 


(4.39) 

(4.40) 

(4.41) 

(4.42) 

(4.43) 

(4.44) 
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4.7  Matrix  Equations  in  the  Traveling  Atomic  Basis 

As  discussed  earlier  if  the  molecular  orbitals  are  expanded  in  terms  of  a basis 
of  static  atomic  orbitals  {xi},  it  would  bring  spurious  coupling  at  large  internuclear 
distances.  To  avoid  this  problem  and  to  account  for  translational  phase  factors,  one 
expands  the  MOs  as  linear  combinations  of  traveling  atomic  orbitals  (TAOs) 

V>7M)  = 5^(r»*)cji(0»  (4.45) 

n 

with  the  TAOs  are  defined  in  chapter  3. 

Operators  can  also  be  expanded  in  a basis  of  TAOs.  The  density  operator  in 
this  basis  is  given  now  in  matrix  notation  by 

P77'  =1  47>P77'<£V  I (4.46) 

and  the  Hamiltonian  matrix  HT)/  is  defined  as 

H77'  = (£7  | H | £7').  (4.47) 

Evaluating  the  matrix  elements  of  f2  and  of  the  electron  kinetic  energy  operator 
K in  the  basis  of  TAOs  and  canceling  terms,  Eq.  4.38  becomes 

iP  = §-1HtP  - PE^S"1,  (4.48) 

where  the  new  overlap  supermatrix  is  formed  by  the  elements 

S’"' = «7  I «V>  = Sirf,  (4.49) 

and  the  Hamiltonian 


H77'  = (KT  + Ve7T)57y+H77; 


(4.50) 
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Now  it  is  convenient  to  define  the  potential  Vqu,  in  Eq.  (2.8),  in  terms  of  an  Ehrenfest 
potential 


Vgu  = VCC(R)  + 


tr[PM] 

tr[P] 


(4,51) 


where  Vcc  is  the  core-core  repulsion  term. 

4.8  Transformation  from  the  Static  to  the  Traveling  Basis  Set 
The  calculation  of  the  matrix  elements  is  first  carried  out  in  the  static  basis 
X and  then  transformed  to  the  TAO  basis  £.  The  transformation  is  done  with  the 
overlap  integrals  between  the  two  bases.  Omitting  provisionally  the  spin  indexes 
7, 7'  and  introducing  a superscript  to  discern  the  x and  £ bases,  one  writes 


P =1  X>P W(*  1=1  «>p<f,«  I (4.52) 


which  yields 

PM  = (S*)-'<x  | OP(£)«  | xXS^r1-  (4.53) 

Here  Sx  is  the  overlap  matrices  given  in  terms  of  the  static  basis. 

4.9  Physical  Properties  of  Interest 
We  are  interested  in  monitoring  the  expectation  values  of  the  electronic  orbital 
(1),  spin(s)  and  total  (j)  angular  momentum  vectors  as  functions  of  time  along 
a particular  trajectory.  To  do  that  let  us  introduce  first  the  expression  for  the 
expectation  value  of  an  general  operator  A,  given  in  terms  of  the  density  matrix  p 

(M*))  = 


tr[p(t)A] 

tv[p(t)] 

tr[E,p  I 6*7 I A] 


377' 


i to  )pv  mu  7'  i 
T,y,p%f(t){to  I A- 1 to) 
Ev'P&WSX 


(4.54) 

(4.55) 


(4.56) 
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The  time  evolution  of  the  Ik  component  ( k — x,  y , z)  of  the  angular  momentum 
is  given  then  by 


Mt))  = 


'EnupH»(t)(&'y'  I h I tui) 


E 


/ i 


(4.57) 


where 


4 = 

/„  = -!«  (zgJ  “ . 

4 = -•'*  (%  - »£)  ■ 


(4.58) 

(4.59) 

(4.60) 


In  the  case  of  a spin  operator,  the  time  evolution  of  the  sx  component  is  given 


by 


(sx(t))  — 


I sx  1 ^7) 


E 


= 2h 


(4.61) 

(4.62) 


for  sy  we  have 


<*»(*)>  = 


I Sy  | f„7> 


Em„  P${t)sZi 


= ±2^ 


i _ E„„  piJ  (t)(&Y  I £m7)(i  -JttO 


(4.63) 

(4.64) 


Eprr'wcn'  ’ 

fJLV  r^V 

where  the  plus  sign  is  for  the  case  7'  = /?  and  7 = ex,  and  the  minus  sign  is  for 
y1  = a and  7 = /3.  For  s2 


&(*)>  - 


1 sz  1 gM7) 

E^ 

_j_  1 ^E/h-  Pflv  (0(£»/  I £m)^77' 
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where  the  plus  sign  is  for  the  case  where  7 = 7'  = a,  and  the  minus  sign  for 
7 = 7'  = p. 

We  are  also  interested  in  calculating  the  total  cross  section  for  the  transition 
j — > j'  at  an  energy  E is  given  by 

1 j j' 

(4-67) 

mj=-j  Tn'.,=-j' 


in  the  above  equation  is  the  fine-structure  cross  section  calculated  with 

and  our  integral  cross  section  <7jmj_>j/m'(  (E)  is 


,,(£0  = 2*  J 


db  b 


,ym\ 


(E,b) 


(4.68) 


for  the  fine-structure  transition  from  a state  jrrij  to  a final  state  j'mly  for  the  pro- 
jectile energy  E,  where  b is  the  impact  parameter  and  Py™*3 is  the  transition 
probability  form  initial  state  jrrij  to  final  state  j'm'j,  at  the  final  time. 

Expressions  for  the  fine-structure  and  total  cross  sections  were  derived  in  Reid.27 
Of  greater  experimental  interest  are  the  multipole  relaxation  and  coherence  transfer 
cross  sections  which  are  discussed  in  detail  by  Wilson  and  Shimoni.82  These  cross 
section  can  be  calculated  from  the  fine-structure  cross  sections  and  the  total  cross 
sections.  For  convenience  we  will  only  summarize  the  expressions  of  those  cross 
sections  which  have  been  measured  experimentally  and  theoretically  calculated  for 
the  Na-He  system.  Using  the  notation  in  Wilson  and  Shimoni,  they  are  given  by 
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When  making  detailed  comparisons  between  theory  and  experiment,  the  cross 
sections  above  must  be  averaged  over  the  velocity  distribution.  The  averaged  cross 
sections  for  a temperature  T are  given  by 


a(T) 


f f(v,T)a(v)dv 
ff(v,T)dv 


where  f(v,T ) is  the  normalized  Maxwellian  energy  distribution  given  by 


f(v,T)  = 4tt 


A4 


3/2 


2i vkBT 


v2ex.p\-iiv'1l(kBT)}, 


(4.74) 


(4.75) 


and  n is  the  reduced  mass  of  the  system. 

4.10  Computational  Aspects 
4.10.1  Spin-Orbit  Integrals 

Integrals  over  the  spin-orbit  (potential)  operator  can  be  factored  into  (vector) 
space  and  spin  parts.  The  computational  choice  to  be  made  is  whether  to  compute 
and  store  integrals  over  molecular  spin-orbitals  or  to  compute  and  store  integrals 
over  molecular  orbitals  of  up  to  three  components  of  the  spatial  part  of  the  operator. 
In  a case  with  no  symmetry,  the  number  of  the  latter  would  be  3/4  the  number  of  the 
former.  With  some  amount  of  spatial  symmetry,  the  ratio  would  be  even  smaller. 
On  the  other  hand,  if  purely  spatial  integrals  are  computed  and  stored,  the  spin 
integrals  of  s must  be  generated  later. 

Our  atomic  functions  are  expanded  in  terms  of  Gaussian  atomic  orbitals  of  the 

form 


g(rx,rix,lx,mx,ax ) = N(nx,  lx,  rnx,  ax)xn^ y1^ z^x  exp  (-axr2x),  (4.76) 

where  the  r*  is  the  position  of  electron  with  respect  to  center  X and  the  normal- 
ization constant  is 

i _ \ {2ax/xfl\AaxYnx+lx+mx)l2 

N(nxJx,mx,ax)  = [(2^  _ 1)!!(2/^  _ 


(4.77) 
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where  we  have  used  the  label  X for  the  centers  A or  B. 
The  spin-orbit  matrix  elements  are  of  the  form 


<9A7  I Hso  | gBl')  = 


VAB 


l l 


E E fr"c 

Z— 1 j m=-l  m'=—l  J 


exp  | lm)(lm  \ 1 | lm')(lm'  \ gB)drc 


(7  | s | 7'),  (4.78) 


where  C can  be  A,  B,  or  a different  center. 

We  describe  below  the  procedure  we  used  to  solve  the  vAB  integral  which  is 
based  on  a procedure  developed  by  Pitzer  and  Winter:62 

Transforming  the  exponentials  of  gA  and  gB  to  center  C in  the  following  manner: 


exp  (- aArA 2)  = exp  ( -aArc 2 - 2a^CA  • rc  - aA  \ CA  |2),  (4.79) 

exp  ( -OLBrB 2)  = exp  ( -aBrc 2 - 2qbCB  • rc  - otB  | CB  |2)  (4.80) 

with 

CA  = C- A,  (4.81) 

CB  = C - B,  (4.82) 


here  A,  B and  C are  the  position  vector  of  centers  A,  B and  C respectively.  Defining 


Dabc  = 47 tN(ua,  lA,  mA,  aA)N{nB , lB,  mB,  aB) 


x exp  (— au  | CA  |2  — aB  \ CB  |2),  (4.83) 
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D 


vab  — 


ABC 

4-7T 


X 


k — — 2(aACA  -t-  a^CB), 

(4.84) 

a = otA  + otB  + /?, 

(4.85) 

kA  = -2c*aCA, 

(4.86) 

kg  — — 2ccgCB. 

(4.87) 

definitions  above  uAB  in  Eq.  4.78  becomes 

a 

Ylf  dr°  \ f d^cxAAylAzAA  ex p(kA 

m=—l  d° 

i-S 
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JS 
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o^ 

1 1 

1 dncxnBBylgz™B  exp  (kB  • rc)Ylm(nc) 

rc  exP  (~arc) 

x (Im  | 1 | lm'),  (4.88) 


expanding  the  exponentials  exp  (k  • rc)  in  spherical  coordinates 

oo  A 

exp  (kA  • rc)  = 47T  ^ Mx(kArc)YXli(9kA,  (pkA)Yx^{6c,  (pc), 

A=0  /i=— A 
oo  A 

exp  (kB  • rc)  = 4tt  Mx(kBrc)YXlt(dkB,  (pkB)YXli(9c,  (pc), 

A=0  /j=— A 


where  Mx  is  a modified  spherical  Bessel  function  of  the  first  kind: 

Mx(x)  = sM ~ 


1 d\X  sinh x 


x dx  I x 


= ixj\{-ix). 


(4.89) 

(4.90) 


(4.91) 

(4.92) 
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Transforming  x A,  Va , za  in  terms  of  CA,  and  separating  variables  of  integrations, 


we  have  for  [a]  in  Eq.  4.88 

= J dQcxAAyAzAA  exp  (kA  • rc)Ylm(nc) 

= f dQc(xc  - CAx)nA(yc  - CAy)lA(zc  - CAz)mA 
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where  we  have  used 


A „ 

E u„( / dacZ‘btfc?cYxAOcAc)Y,m(nc)  (4.93) 
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Calling 
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\ ^ 


/i=-A 

Eq.  4.88  will  become 
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where  the  radial  integral  Q1?-, - is  given  by 


POO 

Qx\  = / drrN  exp  (—ar2)Mx(kAr)Mx(kBr). 
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(4.97) 


82 


Integrals  QN  and  ttabc  are  solved  with  the  procedures  presented  in  appendix  C.  The 
details  of  the  structure  of  the  program  are  presented  in  appendix  A. 

4.10.2  Structure  of  the  Program 

We  start  with  the  pseudopotential  code  that  includes  spin-orbit  coupling.  There, 
we  do  not  need  to  worry  much  about  the  basis  set  given  in  terms  of  coupled  or  un- 
coupled states  since  the  energy  results  are  independent  of  this  choice. 

The  implementation  of  the  Eik/TDMO-PP-SO  (eiktdmoppso)  followed  the 
equations  presented  in  previous  sections  and  what  we  show  next  is  a description 
of  some  details  that  only  appeared  when  implementing  the  code  with  spin-orbit 
coupling. 

In  appendix  A we  present  a diagram  of  the  algorithm  of  the  current  version  of 
the  eiktdmoppso  code  along  with  a detailed  explanation  of  the  subroutines  invoked. 

4.10.3  Choice  of  Initial  Conditions 

For  the  study  of  fine-structure  transitions  we  need  to  start  our  system  in  one 
of  the  fine-structure  states  given  in  terms  of  coupled  basis  | Isjrrij).  But,  as  we 
mentioned  before,  our  code  has  been  written  to  propagate  the  matrices  in  terms 
of  Cartesian  function  basis.  Therefore,  when  we  initialize  our  calculations,  we  first 
create  an  initial  density  matrix  in  the  coupled  basis,  and  we  then  convert  it  from  that 
basis  to  the  Cartesian  basis  (which  is  the  basis  of  the  propagation)  and  finally  when 
the  calculation  ends,  we  follow  the  opposite  by  converting  the  Cartesian  density 
matrix  back  to  the  coupled  one. 

The  procedure  we  created  to  do  the  transformation  of  the  density  matrix  be- 
tween the  basis  sets  is  presented  below. 

The  density  matrix  of  the  system  can  be  expressed  in  terms  of  the  molecular 
orbitals  ip,  the  basis  of  Cartesian  atomic  functions  £c  or  the  basis  of  coupled  atomic 
functions 

V’>p*<V’  i = i e>pc(?c  i=  i (*)&<(* 


P 


(4.98) 
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The  first  step  in  a calculation  is  to  select  the  initial  state  of  the  system  We  do 
that  by  assigning  1 to  the  Pj^  element  of  the  density  matrix.  We  then  transform 
the  density  matrix  from  the  xp  to  £c  basis  that  is  the  one  of  the  dynamics.  Eq.  4.125 
shows  how  it  is  done.  In  what  comes  next  we  explain  how  we  come  up  with  the 
relation  shown  in  4.125. 

Let  us  begin  by  calling  the  Cartesian,  spherical  and  coupled  basis 

respectively.  We  can  write  the  matrix  relations  among  them  as 


ID  = 

| f )TSC 

(4.99) 

ID  = 

| C )Tsj 

(4.100) 

ID  = 

1 Dtjc  = | D(tD-1tsc> 

(4.101) 

where  Ta6  is  the  transformation  matrix  from  coordinates  a to  coordinates  b. 

The  density  operator  p can  be  obtained  in  terms  of  the  £c  and  the  basis  sets 
as 


p = i e )pc((e  i = i e>p’<e  i ■ (4.102) 

Projecting  on  the  left  of  4.102  by  (£c  | and  on  the  right  by  | £c)  we  get 


r 1 e>pj<e  1 n 

(4.103) 

r 1 f )(Ti')-'pi  [(Ti')-1]'  1 r) 

(4.104) 

st(TK)-ipi  [(Ti<)-ijtS* 

(4.105) 

(Tit)-ipi((T,t)-1)t 

(4.106) 

(4.107) 

The  density  operator  p can  be  written  in  terms  of  the  and  the  xp  basis  sets  as 


p = 1 1 = 1 e)pi<e  1 . 


(4.108) 
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Projecting  on  the  left  of  4.108  by  (£J  | and  on  the  right  by  | we  get 

(tj  | | ? ) = SjPjSj  (4.109) 


with  | ip)  = | 


SjCjp'l’Cj]Sj  = SJ~PjSJ 


(4.110) 


CipV’CJ't  — P-? 


(4.111) 


using  in  4.106  the  relation  for  PJ  in  3.56  we  have 

C^pV'c^  = (4.112) 

and  solving  for  Pc  we  get 

pc  = (Tic)-1cJ’P^C,’t((T,'c)-1)t.  (4.113) 

To  obtain  CJ  we  start  with  the  time  independent  Schrodinger  equation 

H | xf>)  =|  ij>) E,  (4.114) 

and  expanding  ip  in  Eq.  4.114  in  terms  of  | ) as  shown  in  3.52,  and  projecting  on 

the  right  by  | we  obtain 

| H | ?)&  = (e'lOC^E  (4.115) 

HjCj  = SjCjE  (4.116) 

multiplying  on  both  sides  of  Eq.  4.116  by  (S-?)-1 

{S’y'H’a  = CjE  (4.117) 

WjCj  = CJ'E  (4.118) 

giving  the  eigenvalue  matrix  equation  in  the  coupled  basis  set.  Therefore  is 

obtained  diagonalizing  W7  in  Eq.  4.118.  The  matrices  and  H7  can  be  obtained 


85 


easily  from  the  corresponding  matrix  calculated  in  the  Cartesian  basis  £c,  to  obtain 


Wj  = (4.119) 

= ((e\cj)r1(e\H\e)  (*-m 

= [({Tjc)-y(Zc  I 0(Tjc)-1]-1  ((T,'c)-1)*(fc  I H I 0(T JC)_1  (4.121) 

= [TJ’c(Sc)-1(T,’c)t]  ((TJC)_1)tHc(TJC)_1  (4.122) 

= TJC(Sc)_1Hc(Tjc)_1  (4.123) 

_ Tjc(Wc)(Tjc)-1.  (4.124) 


Once  the  propagation  in  the  Cartesian  basis  ends,  we  use  the  following  transforma- 
tion to  analyze  the  results, 

P*  = Cj_1TJCPc(Tjc)t(Cjt)~1.  (4.125) 

4.11  Results  and  Discussion 
4.11.1  Pseudopotential  Results 

To  test  the  validity  of  the  method  and  the  accuracy  of  the  computer  code  we 
have  calculated  the  spin-orbit  coupling  splitting  of  Li  and  Na.  Results  are  reported 
in  table  4-2. 


Table  4-2:  Spin-orbit  splitting  for  2P  states  (cm  1). 


Atom 

Hso 

Exp.38 

Li 

0.4 

0.3 

Na 

16 

17 
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We  have  also  calculated  the  diatomic  potential  energy  curve  for  Li-He,  Li-Ne  , 
Na-He  and  Na-Ne.  Figure  4-1  shows  our  results  for  Na-He  and  a comparison  with 
the  ab-initio  results. 

Energies  for  both  the  atomic  and  molecular  energies  were  obtained  with  the 
small  (4s4p/2s2p)  basis.  Despite  the  small  size  of  this  basis  set,  the  results  are  in 
good  agreement  with  the  experimental  and  theoretical  data. 

4.11.2  Pictures  of  the  Evolution  of  Electronic  Angular  Momenta 

To  get  a complete  picture  of  the  collisional  process,  we  examine  the  orbital 
probabilities  as  functions  of  time.  In  figure  4-2  we  show  the  time  evolution  of 
Na(2Pj>mj)  states  for  a particular  trajectory  that  corresponds  to  the  initial  state 
| j,rrij)  = | 3/2, 3/2),  v = 9 x 104  au,  and  an  impact  parameter  of  10  au.  All 
the  populations  oscillate  in  the  time  region  of  20000-45000  au  until  they  come  to  a 
constant  value. 

A simpler  picture  of  the  collision  process  can  be  obtained  by  monitoring  the 
expectation  values  of  the  electronic  angular  momentum  vectors  (1),  (s),  and  (j)  as 
functions  of  time  along  a trajectory.  Each  element  of  these  vectors  is  given  by  the 
expectation  values  in  Eqs.  4.9,4.62-4.66,  where  lx  is  the  matrix  representing  the 
component  of  the  electronic  angular  momentum  operator  along  the  rotating  x axis 
in  the  Born-Oppenheimer  basis.  Using  the  density  matrix  expression  these  matrices 
are  easy  to  derive  and  expectation  values  can  be  computed  at  regular  time  steps 
along  the  trajectory. 

The  exact  evolution  of  the  angular  momentum  vectors  is  complicated,  because 
it  is  determined  by  the  relative  magnitudes  of  three  competing  effects:  the  electro- 
static fields  of  He,  spin-orbit  coupling  in  Na,  and  the  rotation  of  the  internuclear 
coordinates  over  the  course  of  the  collision. 

In  figures  4-3-4-5  we  present  typical  trajectories.  The  He  atom  is  located  at 
the  origin  of  the  space-fixed  coordinate  system.  The  Na  atom  is  traveling  in  the  2 
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Distance  (A) 

Figure  4-1:  Diatomic  potential  energy  curves  of  Na(2P)-He,  including(solid 
line)  and  neglecting  (dotted  line)  spin-orbit  interactions.  Ener- 
gies are  relative  to  the  asymptotic  values  neglecting  spin-orbit 
interactions.  I : Our  results,  II : Ab  initio  results. 


Orbital  Population 
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Figure  4-2:  Population  of  the  Na (2Pj,mj)  states  as  a function  of  time. 
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direction.  We  chose  initial  conditions  so  that  we  could  easily  see  the  time  evolution 
of  (1),  (s),  and  (j).  We  have  used  a moderately  large  impact  parameter  ( b = 10 
au),  v — 9 x 104  au  and  straight  line  trajectory.  The  collision  begins  with  excited 
Na  in  the  2Pz/2,3/2  state.  Our  initial  conditions  do  not  necessarily  correspond  to 
any  easily  attained  experimental  situation.  In  figure  4-3  we  see  the  evolution  of 
(1)  which  seems  to  undergo  precession  about  the  internuclear  axis,  while  during  the 
same  time  as  shown  in  figure  4-4  (s)  rotates  slowly  by  around  90°.  This  results  in 
an  overall  change  in  the  relative  orientation  of  the  two  vectors  and  thus  a change  in 
(j)  as  shown  in  figure  4-5.  Therefore  a fine-structure  transition  has  taken  place. 

In  the  Na*-He  system  jrrij  — »•  j'm'y  transitions  occur  by  the  following  mecha- 
nism. As  shown  in  the  figure  the  system  begins  in  the  j— 3/2  state,  in  which  (1)  and 
(s)  are  aligned  to  each  other.  Then  as  a result  of  the  collision,  the  orientation  of 
both  have  changed.  It  follows  that  both  the  magnitude  and  direction  of  (j)  are  also 
changed  by  the  collision.  Change  in  the  orientation  of  (j)  corresponds  to  rrij  — > rn'- 
transitions.  Usually  a change  of  the  magnitude  of  (j)  represents  a j — > j'  transitions. 
However,  we  must  be  a little  cautious  with  this  last  statement:  actually  it  is  the 
change  of  (j2)  that  corresponds  to  j — > j'\  since  (j2)  (j)2,  it  is  sometimes  possible 

for  the  magnitude  of  (j)  to  change  even  when  (j2)  does  not  change. 


In  tables  4-3  and  4-4  we  present  results  for  the  partial  cross  sections  ►j'm'., 
for  velocities  corresponding  to  two  different  temperatures  400  K and  450  K.  The 
calculations  were  done  at  the  most  probable  velocity  v*  for  the  temperature  T, 
given  by  the  expression 


4.12  Fine-Structure  Cross  Sections 


(4.126) 


where  //  is  the  reduced  mass  of  the  Na-He  system.  These  cross  sections  were  obtained 
by  numerical  integration  over  a grid  of  80  impact  parameters  from  0.0  au  to  20  au 
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Figure  4-3:  Evolution  of  (1)  for  a straight  line  trajectory.  The  initial  state  is 
Na(2P3/2,3/2) 
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Figure  4-4:  Evolution  of  (s)  for  a straight  line  trajectory.  The  initial  state 
is  Na(2P3/2,3/2) 
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z 


Figure  4-5:  Evolution  of  (j)  for  a straight  line  trajectory.  The  initial  state  is 
Na(2P3/2,3/2) 
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Table  4-3:  Calculated  fine-structure  cross  sections  0j,mj-+j'm'.l  at  ^*(400  K) 


j'ytn'j, 

j )mj 

1/2, 1/2  1/2, -1/2  3/2, 3/2  3/2, 1/2  3/2, -1/2  3/2, -3/2 

1/2, 1/2 

668.93 

15.61 

9.84 

13.12 

20.47 

18.81 

1/2.1/2-1/2 

13.37 

636.17 

37.07 

17.30 

37.23 

3.84 

3/2, 3/2 

16.19 

26.44 

635.48 

26.85 

15.50 

25.98 

3/2, 1/2 

17.52 

28.98 

22.59 

632.24 

9.28 

7.66 

3/2, -1/2 

18.40 

17.11 

6.82 

18.21 

660.89 

14.86 

3/2, -3/2 

26.12 

20.23 

15.91 

16.36 

22.68 

641.85 

with  different  step  sizes.  The  tolerance  set  10~7  - 10-5  for  low  and  high  respectively 
seemed  to  converge. 

In  figure  4-6  we  present  the  final  populations  of  some  of  the  Na(2Pjim>)  states 
as  a function  of  the  impact  parameter  for  an  initial  velocity  v*(450K)  = 6.74  x 104 


au. 

We  have  also  calculated  the  total,  multipole  relaxation  and  coherence  cross 
sections.  To  compare  our  results  with  experimental  and  theoretical  results  we  have 
averaged  our  cross  sections  over  a Maxwell  distribution  for  temperatures  of  400  and 
450  K.  In  table  4-5  our  calculated  thermally  averaged  total,  multipole  relaxation  and 
coherence  cross  sections  are  compared  with  experimental  data  and  past  theoretical 
results  including  recent  ones  from  Leo  et  al.73  Our  results  are  similar  in  many  cases 
to  the  past  theoretical  results  but  lie  above  most  of  the  calculations  of  Leo  and  below 


most  measurements. 
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Figure  4-6:  Transition  probabilities  vs.  impact  parameter  in  a Na3/2.3/2+He 
collision,  v * (450ff)  = 6.74  x 104  au 
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Table  4-4:  Calculated  fine-structure  cross  sections  at  v*(450  K) 


j'M, 


3,™*] 

1/2, 1/2 

1/2, -1/2 

3/2, 3/2 

3/2, 1/2 

3/2, -1/2 

3/2, -3/2 

1/2, 1/2 

665.24 

16.16 

17.96 

14.02 

22.06 

24.72 

X/2,1/2-1/2 

16.13 

653.74 

11.58 

35.85 

11.75 

35.85 

3/2, 3/2 

19.22 

39.44 

625.79 

27.70 

19.47 

21.92 

3/2, 1/2 

18.59 

33.94 

20.38 

662.23 

6.40 

4.84 

3/2, -1/2 

28.88 

223.22 

14.96 

16.30 

643.52 

25.27 

3/2, -3/2 

27.66 

15.19 

16.60 

16.84 

17.40 

642.43 
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Table  4-5:  Theoretical  and  experimental  results  for  the  total  and  multipole 
relaxation  cross  sections  (in  A2)  at  400  and  450  K 


a 

T(K) 

Experiment 

Past  theory 

Past  theory73 

Present 

&1/2-+3/2 

450 

106.9±1183 

75. 227 

85.4 

77.55 

400 

86  ±884 

00 

CO 

OO 

81.5 

77.05 

cr3/2->l/2 

450 

46±484 

43. 627 

44 

46.46 

400 

57±583 

43. 838 

44.5 

46.39 

adep 

450 

13.2±1.883 

12.5585 

10.8 

16.07 

400 

15.6±1.686 

10.982 

11.0 

16.12 

o\ 

450 

133±783 

119.785 

103.1 

108.87 

2 

400 

90±986 

11282 

103.5 

108.45 

cr\ 

450 

112±683 

110.085 

88.6 

90.84 

2 

400 

128±1387 

9882 

90.2 

90.89 

_co/i  tr 
u 1 

450 

-66.3±10.383 

-7585 

-57.5 

-60.39 

2 

400 

-64±  6.486 

-62. 482 

-56.5 

-59.83 

_co/i  tr 
V 3 

450 

-6.5±1.183 

-6.985 

-5.77 

-10.36 

2 

400 

-10.15 

CHAPTER  5 
CONCLUSION 

The  purpose  of  this  study  has  been  to  develop  a time-dependent  methodology 
to  explore  in  detail  the  dynamics  of  alkali  atoms  colliding  with  rare-gas  atoms  at 
two  very  different  energy  regimes: 

• A high  energy  (0.5  to  30  keV)  regime  to  study  the  time-dependence  of  the 
electronic  excitation  of  the  alkali  atom. 

• A low  energy  (hyperthermal)  regime  with  energies  of  the  order  of  fractions  of 
one  eV  to  study  the  spin-orbit  recoupling  of  an  initially  excited  alkali  atom. 

A new  formalism  which  combines  pseudopotentials,  spin-orbit  coupling  potentials 
and  the  eikonal  time-dependent  molecular  orbital  method  has  been  developed  to 
study  the  dynamics  of  the  active  electron  and  nuclei  by  reducing  our  many  electron 
systems  into  one  active  electron  systems.  The  conclusions  of  this  work  are  given  in 
the  same  order  as  that  of  the  preceeding  chapters. 

5.1  Pseudopotentials 

An  overview  of  the  pseudopotential  theory  has  been  presented.  The  approach 
we  chose  in  this  work  was  the  /-dependent  pseudopotential  method  including  core- 
polarization suggested  by  Preuss  and  collaborators.58  This  choice  was  based  on 
the  following  premises:  a)  it  reproduces  accurately  the  interatomic  potentials  of 
the  molecular  systems  of  interest;  b)  all  the  pseudopotential  parameters  and  basis 
sets  were  available  in  the  literature;  c)  the  pseudopotential  operators  are  expanded 
in  terms  of  standard  functions,  for  which  efficient  analytical  methods  to  solve  the 
integrals  are  available  in  the  literature. 

We  successfully  reproduced  the  experimental  ionization  energies  of  the  Li  and 
Na  atoms.  We  calculated  the  potential  energy  curves  of  the  Li-He,  Li-Ne,  Na-He 
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and  Na-Ne  diatomic  systems  and  found  them  to  compare  very  well  with  electronic 
structure  calculations. 3’4,31,41’67  We  have  them  available  for  a series  of  alkali  and 
rare-gas  atoms. 

We  found  that  the  basis  set  size  has  an  important  impact  on  the  calculations 
and  that  convergence  with  respect  to  potential  energies  versus  distances  is  reached 
only  when  basis  sets  containing  several  d- functions  are  used. 

We  conclude  from  our  results  that  /-dependent  atomic  pseudopotentials  are 

very  convenient  and  accurate  for  studies  of  an  alkali  atom  with  a rare-gas  atom. 

Our  present  results  showed  that  they  give  excellent  potential  curves,  by  comparison 

with  previous  work,  and  also  led  to  the  correct  non-adiabatic  couplings  between 

electronic  states.  The  potentials  we  obtain  are  of  sufficient  quality  so  that  they 

could  also  be  used  to  describe  the  spectra  of  alkali  atoms  in  a rare-gas  environment 

5.2  The  Eikonal-Time  Dependent  Molecular  Orbital  Method  Including 

Pseudopotentials 

A novel  treatment  based  on  the  Eik/TDMO  method  in  terms  of  pseudopoten- 
tials was  introduced  to  study  the  time  evolution  of  electronic  excitation  in  collisions 
of  alkali  atoms  with  rare-gas  atoms,  and  was  called  the  Eik/TDMO-PP  method. 

The  calculation  of  time-dependent  properties  such  electronic  populations  and 
polarization  parameters  for  these  systems  provided  a unique  insight  into  the  elec- 
tronic rearrangement  mechanism.  It  was  noted  that  electronic  populations  tended 
to  fluctuate  during  collisions  and  to  settle  at  a final  value.  We  found  that  the  final 
value  of  the  electronic  populations  was  a function  of  the  projectile  energy  and  the 
impact  parameter. 

We  calculated  the  integral  cross  sections  for  the  electronic  excitation  from  the 
ground  state  of  Li  and  Na  to  their  first  excited  state  in  collisions  with  He  and  Ne  in 
the  range  from  0.5  to  30  keV  depending  upon  the  system.  We  obtained  reasonable 
agreement  with  measurements  for  those  systems.  From  our  results  we  conclude  that 


99 


several  atomic  d-orbitals  were  needed  to  describe  2s  to  2p  transitions  of  Li  and  the  3s 
to  3p  transitions  of  Na  and  that  a large  basis  set  was  necessary  for  obtaining  accurate 
integral  cross  sections.  We  found  that  the  basis  set  size  also  affected  probabilities 
versus  impact  parameters  and  therefore  the  differential  cross  sections. 

5.3  Spin- Orbit  Recoupling 

We  extended  the  Eik/TDMO-PP  formalism  to  include  spin-orbit  coupling  for 
the  study  of  angular  momenta  recoupling  during  hyperthermal  alkali-rare-gas  atom 
collisions.  The  work  involved  the  use  of  pseudopotentials  and  it  also  required  the 
construction  of  the  molecular  spin-coupled  states  from  atomic  spin-coupled  states. 

The  calculation  of  the  spin-orbit  Hamiltonian  was  possible  using  a modified 
version  of  a procedure  available  in  the  literature.62  We  calculated  the  spin-orbit 
coupling  splittings  for  Li  and  Na,  and  our  results  compared  well  with  experimental 
data.  We  generated  the  potential  energy  curves  of  Na-He  including  spin-orbit  and 
they  were  in  good  agreement  with  other  theoretical  results. 

We  presented  pictures  of  the  evolution  of  the  expectation  values  of  angular 
momentum  vectors  (1),  (s)  and  (j).  They  allowed  us  to  interpret  the  results  and  to 
obtain  information  about  the  mechanism  leading  to  the  fine-structure  transitions. 
We  found  that  the  evolution  of  (1)  underwent  precession  about  the  internuclear  axis 
while  (s)  rotated  slowly.  The  change  of  orientation  of  (1)  relative  to  (s)  resulted  in 
a change  of  (j),  i.e.,  a fine-structure  transition. 

Finally,  we  calculated  the  fine-structure  cross  sections  of  Na-He  at  400  and 
450  K.  Our  calculated  total,  multipole  relaxation  and  coherence  cross  section  were 
compared  with  experimental  and  theoretical  data.  Our  results  were  similar  in  many 
cases  to  theoretical  results  and  showed  correct  magnitudes  and  trends  compared  to 
experiment.  Theoretical  results  are  consistent  among  themselves  but  generally  they 
were  below  most  measurements.  The  agreement  in  general  was  acceptable. 
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5.4  Computational  Aspects 

We  coded  a new  pseudopotential  program  entitled  eiktdmoppso  entirely  in 
Fortran  90.  We  presented  a complete  derivation  of  the  hamiltonian  integrals  in  the 
chapter  2 and  in  the  appendices  A,B  and  C.  We  rewrote  in  Fortran  90  The  eiktdmo 
to  efficiently  calculate  all  the  necessary  components  appearing  in  the  integration  of 
the  equations  of  the  time-dependent  density  matrix  formalism  including  pseudopo- 
tential and  spin-orbit  coupling.  We  included  in  appendix  A a flow  diagram  of  the 
time-dependent  density  matrix  computer  code. 

It  is  remarkable  that  our  theoretical  and  computational  treatment  of  the  elec- 
tronic and  nuclear  dynamics  can  provide  the  information  needed  for  optical  transi- 
tion phenomena  at  both  hyperthermal  and  high  collision  energies. 

5.5  Future  Work 

The  formalism  and  computer  implementation  of  the  method  are  suitable  for 
completing  studies  similar  to  the  ones  presented  above,  on  other  alkali-rare-gas  atom 
pairs,  as  well  as  the  alkaline-earth-rare-gas  atom  pairs.  The  present  work  will  allow 
us  to  study  spectra  and  dynamics  in  clusters  containing  an  alkali  atom  or  alkaline- 
earth  ions  and  several  rare-gas  atoms. 


APPENDIX  A 

FLOW  CHART  FOR  THE  EIKTDMOPPSO  PROGRAM 
A.l  Program  Features 

The  program  eiktdmoppso  is  the  latest  extension  of  the  original  eiktdmo  pro- 
gram (written  by  Keith  Runge70).  Most  of  the  subroutines  of  the  original  code  have 
been  upgraded  into  fortran  90,  some  have  been  rewritten  and  many  more  have  been 
added.  As  a result  the  eiktdmoppso  code  contains  several  new  options. 

The  type  of  calculation  the  program  will  perform  is  controlled  by  a single  option 
number  (read  from  the  input)  called  ’task’.  These  task  numbers  go  from  [1-15];  task 
numbers  [1-9]  tell  the  program  to  run  the  dynamics  of  the  system  using  a particular 
potential  choice,  as  given  below  (task  numbers  [1-5]  are  the  options  from  the  old 
eiktdmo  program): 

1.  Straight  lines  (V^u=0). 

2.  Coulomb  trajectory  (no  electronic  contribution). 

3.  Screened  coulomb  potential. 

4.  Average  effective  potential. 

5.  Ehrenfest  potential. 

6.  Pseudopotential  calculation. 

7.  Spin-orbit  coupling  calculation. 

8.  Pseudopotential  calculation  using  parameterized  hamiltonian. 

9.  No  available  yet. 

Task  numbers  [10-15]  tell  the  program  to  diagonalize  the  W matrix  and  depending 
on  the  number  it  will: 

10.  Print  eigenvectors  and  eigenvectors  of  the  W matrix. 

11.  Print  electronic  energies  neglecting  core-core  repulsion  contribution. 
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12.  Print  electronic  energies  including  core-core  repulsion  contribution. 

13.  Print  electronic  energies  including  core-core  repulsion  contribution  and  spin- 
orbit  coupling. 

14.  Print  electronic  energies  using  parameterized  Hamiltonian 

15.  Print  eigenvalues  and  eigenvectors  from  calculation  including  spin-orbit  cou- 
pling. 

A. 2 Subroutine  Description 

In  this  section  we  turn  our  attention  to  the  description  of  the  features  of  the 
main  subroutines. 

• read_input  .f  90:  Read  the  input  in  either  the  old  (eiktdmo)  style  or  the  new 
(eikdtmoppso)  one.  Calculate  factors  to  be  used  throughout  the  calculation 
and  also  the  normalization  constants  of  the  gaussian  primitives.  If  requested, 
print  out  the  information  read  from  the  input. 

• get_overlap.f90:  Calculate  the  overlap  matrix  S at  t0  as  well  as  the  matrices 
S1/2  and  S-1  which  are  needed  for  the  transformation  between  basis  sets. 

• init  .f  90:  Allocate  matrices  and  initialize  them  as  well  as  the  other  variables. 

• print_energies . f 90:  (Invoked  only  when  task  > 10)  Print  electronic  energies 
as  explained  above. 

• set_p0.f90:  At  t0  construct  the  density  matrix  in  terms  of  the  basis  set  and 
transform  it  to  the  cartesian  set. 

• set-t0. f 90:  Save  the  information  of  the  variables  at  to  in  case  the  current 
step  is  rejected. 

• buildq.f90:  Compute  the  change  in  the  density  matrix  Q. 

• get_p0.f90:  Calculate  the  current  value  of  the  reference  density  PO. 

• test_step . f 90:  Compute  the  quotient  ||  Q ||  / ||  PO  ||  and  compare  it  to  the 
tolerances.  Then  decide  whether  to  accept  or  reject  the  step. 

• update. f 90:  Update  the  variables  and  matrices  if  step  is  accepted. 
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• pot. f 90:  Calculate  the  current  potential. 

• props. f 90:  If  requested,  calculate  and  print  time-dependent  properties  (e.g. 
Lowdin  populations,  alignment  parameters,  etc.). 

• out. f 90:  Report  final  results. 

A. 3 eiktdmoppso  Flow  Diagram 

We  present  now  the  flow  diagram  of  the  program;  it  has  been  split  in  three 
parts  for  convenience. 
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Figure  A-l:  Cycle  for  printing  energies  and  initializing  dynamics. 
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Figure  A-2:  Relax-and-drive  cycle. 
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Figure  A-3:  Update  and  continue. 


APPENDIX  B 
SUBROUTINE  POLAR. F90 


B.l  Subroutine  polar. f 90 

This  program  evaluates  the  matrix  elements  of  the  polarization  Vpoi  and  cross 
terms  Vec  of  the  potential  energy  operators  in  terms  of  cartesian  gaussian  functions 
using  a procedure  presented  by  Schwerdtfeger  and  Silberbach  .63 

The  procedure  has  the  advantage  of  reducing  the  integral  to  standard  functions, 
including  the  error  function  erf(x),  the  Dawson  function  d(x),  and  the  Dawson-error 
hybrid  function  m(x): 


B.1.1  Evaluation  of  the  rK  Integral 

Without  loss  of  generality,  we  want  to  solve  the  integral 

A02  = ( 9i  | xmi ym2 z™3 r~k (1  - \ g2)  (B.2) 


with  CGTO’s  gi  and  g2  at  centers  1 and  2: 

9j  = (x-  Xj)^(y  ~ VjT 23 (z  - exp  (— /?j(r  - r,)2),  j = 1, 2.  (B.3) 

Here  k > 2,  q G N,  e ff0,  a,  ^ G K+,  the  angular  indexes  i — 1,2,3,  and  the 
center  indexes  j = 1,2.  The  three  different  centers  are  also  denoted  in  the  indices  of 
7102  with  0 representing  the  origin.  For  the  following  we  omit  the  subscripts  of  I. 
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Let  us  start  by  applying  the  binomial  formula  to  the  cutoff  function,  as  well  as 
to  the  x,  y and  z factors  of  the  CGTO’s: 

9 


(1  - e~ar)q  = ( 9 ) (— 1 Ye-»ar*  (B.4) 

(x  - Xj)"j  = (-1  )"*-*« xJtf-War'W.  (B.5) 

tki=o  Wi' 


We  can  rewrite  Eq.  B.2  as 

1 = 'E,  f in,  n) 

9 

x f xllylhl*  exp[-/?i(r  - rx)2  - /?2(r  - r2)2]  V (-l^e^’dr. 

J® 3 t^o  W 

Note  that  the  first  summation  is  taken  over  six  indices  rjij  e {0, 1, ... , ntj } : 

nn  «2i  «32 


and  / is  defined  as 


E-E  E-E 

9 911=0  921=0  932  =0 


1 ,j  v*i/ 


(B-6) 


(B.7) 


where 


(£ij  > %2j  > £3 j ) — (xj > Uj  1 zj ) (B.8) 

and  l1-  — r/n  + rji2  + mj  G N0.  Defining  the  Laplace  transformation  to  remove  the 
denominator  r~k(k  G N): 


r-1(fc/2)  [ e-tr2tk/2~ldt, 

J R_|_ 


(B.9) 
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where  T is  the  gamma  function.  Using  this  transformation  the  above  integral  can 
be  written  as 


I=T~'(k/2)Y,fM  E 


V R3xR+  A*=0 

x exp[— /3i(r  — iq)2  — Pi(r  — r2)2  — (ya  + t)r2]dtdr.  (B.10) 


Using  the  lemma  of  Fubini  for  infinite  intervals 


/ 

Ju3 


dtdr 


I 

J ]R-j 


■dr  dt  , 


(B.ll) 


'R3xR+ 

we  can  first  take  the  integration  over  the  range  r G M3 . Before  doing  so,  we  transform 
the  exponent  in  B.10 

n2  b2 


/5i(r  — ri)2+/?2(r  — r2)2 + {ya  + t)r2  = 
using  the  shorthand  notations 


V^r 


yjd’tpi . Qf/i 


+Pirl+(32r2,  (B.12) 


+ t — (3 i + /?2  + /-tor  + t,  b — (&i,  b2,  63)t  — /?iri  + /?2r2,  b — \ b | . (B.13) 
By  formal  rearrangement  we  get 

I = r-\k/2)e-^-^  V /(Ty,  n)  f Y]  (q)  (-1  )^/2-ig62/o^ 


/i=0 

x / xl'yl3zl3  exp 
I R3 


/ 

JR3 


b 1 

2' 

drdt.  (B.14) 


We  try  now  to  solve  the  integral  [a]  on  r in  Eq. 

/ 

./R3 


/•»  />> 

x ly  2z3  exp 


1 n 2 


- 


\Zatn . . 


dr. 


(B.15) 
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Let  us  split  the  last  integral  in  its  x,  y,  z components.  Substituting  u>  = y/a^r  — 
b / y/atil,  if  we  take  the  x component  we  get  (calling  it 

\ i? 


a,x 


a,x 


= —[( 

y/at)i  J \ \fatn 


U>X  + bx  ) e ^Xdujx 


(B.16) 


calling  ujx  = u,  bx  — b,  Z"  = Z,  yfa^  = a we  have: 

= - [ (~  + b)1  e-“2duj. 
a J \a  / 


(B.17) 


Using  the  binomial  theorem 


a,x 


(B.18) 


being  lj  an  axis,  when  i is  odd  the  integral  goes  to  zero,  so  for  even  values  of  i we 
get 


a,x 


-iSCMi)1/-1-'* 


(B.19) 


therefore 


w ^ 


/^V  (» - 1)!! 

2*/2 


^=^E  (») 


a * — ' \t 

i=0 


(B.20) 


valid  for  Z"  ^ 0.  Using  the  convention  (—1)!!  = 1,  this  leads  to 


I = Y.  /(I. »)  E 9(C.  n 

n C 

x / (B.21) 

"'K+  /j=0 

The  second  summation  is  taken  over  the  three  indices  Q G {0, . . . , Z?},  i = 1, 2, 3 : 

J?  ^ 

= £ (Ci,  C2,  Co  even),  (B.22) 

C Ci=o  C2=o  Cs=o 


Ill 


and  g is  defined  as 


»^)=n(3^ 


(B.23) 


The  are  defined  in  Eq.  B.1.1,  and  + 1%  + I3  ~ (Ci  + C2  + Ca)/2  € No-  Now 

we  turn  to  the  last  integral  over  t G K+.  As  we  will  see  below,  the  integral  and 
sum  (over  fi)  are  interchangeable  only  for  the  conditions  mentioned  before.  For  the 
general  case  (2s  - k G Z)  we  have  to  solve  the  integral 

J = [ vf9)(-l Y{a^t)-s-z'2tkl2-1eb2l^+t)dt 

JR+  p \Av 

= limV(^V-l )•*  [ (a„  + 1)— s/2^/2-1c(*,/a^+*)<ft>  (B.24) 

£HOO/^  W •/[<>,«] 

where  s = sjl  G N0,  a^,  b G K+.  We  must  consider  two  cases  whether  the  power  of  A; 
of  rk  is  even  or  odd. 

Case  1.  k = 2j,  j G N.  We  substitute  = b2/(a^  + t),  and  make  use  of  the 
binomial  theorem  for  the  factor  t: 

j- 1 


. 3 i2\  b2 

xdr)[s-j  + v+-,—  ),r]= — (B.25) 


2 a, 


a,,  + e 


with 


dv(a,x)  = j wa  le“duj.  (B.26) 

J [77, x] 

For  positive  first  arguments  of  dv,  the  limit  and  sum  in  Eq.  B.25  are  interchangeable 
and  the  d function  may  also  be  written  as  an  incomplete  7 function  [d(a,  x)  — d0(a , x) 
if  a > 0] 

7 (a,x)  — f e_tta_1df  = (— l)ad(a,  — x)  (B.27) 

J[  0,*] 
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with  a E R+.  From  a p-fold  integration  by  parts  we  get  the  general  formula  ( p 
cl  — 1/2  E No ) 


(2p-l)H 


p-i 


d(v  + 1/2, x)  = k ^V(-l)pe*  ^(v^)  ~ EGl)1^1'2^!)!!  I ' (B28) 

where  d(x)  is  the  Dawson  function  introduced  in  Eq.  B.l.  By  convention  the  sum 
in  B.28  vanishes  if  p = 0.  To  find  a formula  for  the  d function  with  negative  first 
arguments  we  set  p — j — s — v — 1 and  integrate  by  parts  p times  using  the  equation 

t~p+ Vv  IX 


j 

J \n 


t-P-'^e'dt  = 

M ~P  + 1/2 


~P+  I/2  /[r,,! 


t-p+l/2etdt 


(B.29) 


For  the  first  term  on  the  right  hand  side  in  B.29  we  have  to  evaluate  the  limit  77  — 0, 
where  77  = 62/(ap  + e)(e  — > 00): 


1 /_\  / b2  \ -p+1/2 


lim  Sq(x)  = lim  V'  (-l)14  ( - 


+ e 


(B.30) 


Substituting  ax  = e + ap  — //a  leads  to 


a \ p— 1/2 

.62 


= (£)'  J™  E +^r1/2- 


P=0 


(B.31) 


Expression  B.31  tends  to  zero  with  x approaching  infinity  only  if  q > p — 1/2  or 
q > &/2  — s — 1 for  all  u.  From  this  we  derive  the  formula 

p 


d(l/2  - p,  x)  = e1  ^2Cppd(^)  - Y CiPxi~P~1'2 
with  coefficients 

Cip  = ri(P+l/2-,)-=2‘(H^f. 


(B.32) 


i/=i 


(B.33) 
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Case  2.  k — 2 j + 1 ,j  £ N.  We  rearrange  B.24  and  substitute  z — b2 /a^ 
62/(a/J  + 1): 

9 /_\  . -,*/■>  />  /b2  b2  \ s+3/2 


J2  = 


lim'T(V)(-l  )-(^)'+3/2/  (*-J—  ) 

' ^[0,77]  a„+V 

fj=0  vrv 

x f (- zj~1^2e~zdz,  (B.34) 

J [o,(V Me]  \an  J 

where  e comes  from  the  left  to  1.  We  must  consider  two  cases,  whether  s — j is 
positive  (including  zero)  or  not. 

Subcase  1.  s-jeN0.  We  make  use  of  the  binomial  theorem,  set  e = 1 (which 
is  allowed  in  this  case)  and  write 

*=E(!W^1(£)^/h^ 

/j=0 

<*»> 


i/=0 


The  incomplete  7 function  was  defined  previously  and  we  can  use  B.27  and  B.28  to 
derive  the  7 function  from  the  erf(x ) B.l: 

p- 1 


7(P  + *)  = ^2P2p_i1^‘-  ( erf(y/x)  - ex  ^ ^+1/2 

V P=0 


2^ 


(2/x  + l)!! 


(B.36) 


with  p £ No.  By  convention,  the  second  sum  in  B.36  vanishes  if  p = 0. 
Subcase  2.  s — j G Z ^ 0.  We  have  to  solve  the  integral 

p-j- 1/2  /.  ,2 

/ zj~l/2(b2 /a  ^ - z)~pe_zdz,  e < — , 

J [0,e] 


5 O' 


/ 

</  [0 


a; 


2j 


,e]  (1  - S2)P 


e^l^^dx  = he(j,p,b2/afi),e  < 1 


(B.37) 
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with  p = j — s.  As  a first  step,  we  extend  the  numerator  in  B.37 


, , . , f X2j  - X2(j  + X2(j  _ 2 

K(3’p’a)  = L — — e 

= he(j  - 1 ,p,a)  - he(j  - l,p  - l,o),  j e N,  p e Z. 


(B.38) 


This  is  a recursive  definition  of  the  general  ht  function  and  we  get  the  he  function  of 
the  zeroth  order,  which  we  denote  by  Ht(p,a)  = he(0,p,a)  by  multiple  application 
of  B.38 

he(j,P,  a)  = (“0  (-1  VHe{p  - rj,  a).  (B.39) 

Note  that  this  integral  exists  only  for  e < 1 if  p < 0 and  we  cannot  exchange 
the  limit  and  the  sum  in  B.34  recursively  in  the  limit  e = 1 if  q > j — 2: 

H(p,  <0  = --^ta>1i3g(P  - 1.  «>  - - 2, 1), 

2(p  — 1)  p - 1 

H(0,a)  = a~1^2erf(y/a), 

H(l,  a)  = 2e~am(y/a).  (B.40) 


We  can  either  use  B.40  for  the  defining  the  function  H or  we  must  write 

,)(-l)"if«(p,a).  (B.41) 

/V 

The  function  m(x)  is  defined  in  B.l.  p — p can  become  negative,  and  therefore  we 
define  the  H function  also  for  negative  first  arguments  similar  to  the  d function  in 
B.l: 


£ 

/i=0 


l(— l)"if(p,a)  = limj^  ( 
£_>  n= o ' 


H(—p,  a)  = f (1  - x2)pe~ax2dx 
h o,i] 

= 5E(^)(-ira-‘'-1/27(‘'+l/2.a). 

Z i/=0  ' ' 


(B.42) 
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We  can  now  write  down  the  solution  for  negative  values  of  (s  — j ): 

j2  = 2±  £ (3)(-l )'H(j  - s - n,b2/a,).  (B.43) 

/j=0  W r,=0  V" 

Now,  the  integral  is  completely  solved  and  we  summarize  the  results  using  Eqs. 
B.21,  B.25,  B.35,  B.43: 


/ = r-1  £/(>>, «)  n Y 

r)  C /*=° 


(B.44) 


with  k = 2ki+k2,  fci  G N,  k2  G {0, 1},  <7  G N,  q > «(&)— X^i=i  K(raj)  — 1,  «(n)  = n/2 
if  n is  even,  and  «(n)  = (n  + l)/2  if  n is  odd,  + (32  + ^a,  b = | fail  + /?2r2  |, 

the  functions 


= b 2i+2j  3 XI 

t/=0 

«i(i.  «> &)  = aJ~s~1eb2/a  X 

i/=0 


O'  - 1 


a \ j+H-1/2 


s — J 


(-ir  ( 


52/ 


7 


1 621 

J + " + 9.— 

2 a 


(B.45) 

(B.46) 


for  s — j G N0 , and 


= 20^— (*)  (-1)"H0  -s-v,^)  (B.47) 

i/=0  ' ' 

for  j — 2 G N.  The  T function  is  defined  in  B.27,  the  H function  in  B.40  and 
B.42,  the  erf,  d,  and  m functions  in  B.l,  and  the  / and  g function  and  appropriate 
sum-conventions  in  B.7  and  B.23. 

B.l. 2 Special  Case  b=0 

This  case  occurs  when  both  CGTO’s  and  the  operator  are  placed  at  the  origin: 


-fooo  — hm/so.B.2) 
&->o 


(B.48) 


which  is  equivalent  to  solving  the  integral 


/ 

J R3 


XPlyP2zP3r  k 


± (:)(-» 
7»=0  Kr/ 


He~Uh+p2+Ha)r2  fa 


(B.49) 
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with  pi  = Hu  + ni2  + rrii,  i G {1,2,3}. 


Case  1 . k = 2 j,  6 = 0.  The  sums  over  ry  and  ( in  B.21  break  down  to  just  one 
term  and  the  integral  is  reduced  to  an  elementary  one.  The  result  is 


/ooo  = r-'GV72  (®) (-!)'•<— 3/2  i>i) 


n=0 


v=0 


v (B.50) 


Using  the  formula 


t H)" 


f- 


^v\(x  + v)  x(x  + l )j 


, x G K, 


where  (—z) u,  z G R,  and  u G N0  is  the  factorial  function,  we  get 

/.»» = ^ ff<»  - ■ + })"'  £ ( !)  (-i  ra-i-w. 

t=0 


2'  ' \p, 

fi=0 


(B.51) 


(B.52) 


Case  2.  k = 2j  + 1,  b = 0.  Using  the  mixed  Gaussian-Laplace  transformation 


r-V-'  = ^=T-\j)  [ e-0i+*?)r (B.53) 

V71-  ./IR*. 

instead  of  B.9,  the  first  integration  leads  to  Eq.  B.50  substituting  only  — > a^+tl- 

We  now  integrate  B.50  over  t2  G R+.  Two  cases  have  to  be  considered,  whether  or 
not  j — s < 0.  In  the  first  case  we  use  the  expression 


-p- 3/2 


2pp! 

ap+1(2p  4-  1)!! ’ p~  S~J 


(B.54) 


and  get 


-fooo  — 


2*(s-j)! 

(2s + 1)!! 


(— l)po 


j-s-l 

P 


(B.55) 


APPENDIX  C 

SUBROUTINE  SHORT  _RANGE.F90 

C.l  Subroutine  short_range.f90 

This  program  has  been  written  to  evaluate  the  matrix  elements  of  the  short- 
range  operator  in  terms  of  cartesian  gaussian  functions61  for  any  l in  either  the 
operator  or  the  gaussian  functions. 

The  semi-local  (/-dependent)  short-range  operator  has  the  form 

V£'(?c)  = X B',< eM-PiA)Pi,c,  (C.l) 

where  Biti  and  Piti  are  pseudopotential  parameters  adjusted  to  experimental  data 
and  Pitc  stands  for  the  projection  operator  on  angular  symmetry  l 


P,,c  = '£\Ylm(rc))(Ylm(rc)\.  (C.2) 

771 

Here  | Yim{fc))  represents  a real  orthonormal  spherical  harmonic  function  cen- 
tered on  the  core  C ( fc  = rc/rc)- 

The  integral  between  gaussian  functions  g a,  and  the  short  range  can  be 

written  as 


7ab  = [ dr°  [ dtlcgAYimi^c)  r2c  exp  f dClcYim(flc)gB  ■ (C.3) 

m=-r 0 U i U . 

The  derivation  of  the  solution  of  the  above  equation  is  similar  to  the  one  pre- 
sented in  the  chapter  on  spin-orbit  coupling.  After  following  the  steps  4.78-4.96,  we 
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get 


TlA  l A ™A  nB  Ib  BIb 

jab  = 4 nDABc  Y1 

a=0  6=0  c=0  d=0  e=0  /= 0 


'nA\  ( Ia\  ( mA\  ( nB\  (Ib\  ( mB' 


)\d 


f 


x CAnxA ~aC  A1  A ~bC  A™  A ~cCBnxB ~ dCBlyB ~eCB?B ~f 


oo  oo 


l 


EE« 


a+b+c+d+e+f+2. 


x>  > 2(kA,kB,a)J2^m^fiL  (C-4) 

A=0  A=0 


m=—l 


the  angular  integral  is  given  by 


n$£,  = E yv(«»)  / ^iaj/li°n„(n)v;m(n), 


(C.5) 


/*=-A 


and  the  radial  integral  is  given  by 


:AA 

r»oo 


roo 

Q\\~  / drrN  exp(-ar2)Mx(kAr)Mx(kBr). 
Jo 

The  only  nonzero  terms  in  the  sum  over  A are 


(C.6) 


max{l  — a — b — c,0)<\<l  + a + b + c (C.7) 

and  likewise  for  A.  Also  l + a + b + c — A must  be  even. 

C.2  Special  Cases 

C.2.1  kA,kB  = 0 

The  exponential  terms  for  both  A and  B is  reduced  to 


exp(O.rc)  = 1 


(C.8) 
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and  the  short  range  integral  is  reduced  to 

i 


Iab 


Dabc 

47T 


r c 


drc 


x r> 


4BC  . 


f dncxn/y^z^Y,m(ac) 

exp(-crj)  17  dncx’'/V,B’z’S’Ylm(nc) 

l 

^ ' ^|M+bl+”MQ”B+!B+’T»B 
m=—l 

The  radial  integral  takes  the  form 

roo 

<5^(0:)=  / drr^exp  (— ar2), 

Jo 

and  the  angular  integral 

= / 

= ±r;zj  M*‘+ryt 


',6+s  ~c+t 


(C.9) 


(C.10) 


(C.ll) 


r,s,t 


C.2.2  or  ks  = 0 

Let’s  choose  ks  = 0 for  the  moment,  the  short  range  integral  becomes 
A B _ A4BC 


47T 


S f drc\  [ dncxnAAyl^z^A  exp(kA-Tc)Yim(nc) 
n=-r  0 ^ 

xr^exp(-arc)  J dDcxnBByl§  z^YimiDc) 


• (C.12) 


Expanding  the  exponential  exp  (kA  - re)  in  spherical  coordinates  and  trans- 
forming xa,Va,  za  to  point  C and  separating  variables  of  integration  we  obtain 

n A Ia  iriA 


Iab  - Dabc  ^ 

a=0  6=0  c=0 


nA\  [Ia  \ ( f^A 
a 


C AnA~aC AlA~bC  AmA~c 

x y z 


b J \ c 

00 

x E£3fWc'M+'+/+2(*U,a)  Y,  (C13) 


A=0  A=0 


m=—l 
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and  the  radial  integral  Qx  is  given  by 


drrN  exp  (—ar2)Mx(kAr). 


(C.14) 


C.3  Evaluation  of  the  Integrals 
C.3.1  Angular  Integrals 

Expanding  the  real  orthonormal  spherical  polynomials  Yx h in  terms  of  x,y  and 

z: 


YXfi  = ^2  y^txrysz\  r + s + t = X 

r,s,t 


(C.15) 


the  complete  angular  integral  will  be 


«£.  = E 

/i=-A 


,r,s,t 


\ l 


r,s,t  u,v,w 

The  evaluation  of  the  integral  is  straightforward 


£ E S/Jz.  I dttxa+r+uyb+s+vzc+t+w.  (C.16) 


47t  f d,Vtxly^zk  = 0,  i,j,k  odd, 

_ (* -l)!!(j -!)!!(* -1)H 


(i+  j + k + 1)!! 


i,j,k  even. 


(C.17) 


C.3. 2 Radial  Integral 
Single  Power  Series 

We  start  with 


POO 

Qxx  = / drrN  exp(-ar2)Mx(kAr)M-x(kBr), 

Jo 

and  doing  a change  of  variables  r = r/y/a  and  dr  = dr/y/a  we  get 

Qxx  = f dr("^)JVexP(_Q!r2)MA^)M^(^)- 

y Ot  J o yd  y Oi  y Ot 


(C.18) 


(C.19) 
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Replacing  only  M\(kA.r / <Jot)  by  its  expansion  in  a power  series  we  get 


1 r°° 

(?X  = -n± t drrNexp(-r2) 

a 2 Jo 


( 

[ kAr 

A oo 

(&) 

Lv^. 

^j!(2A  + 1 + 2j)H  J 

k A «o  { ~ (^y  ^ 

= ^exp(-^)  E,-!(2AV !;«)!! 


Mj(tBr) 


k. 


A °° 


E 


/kiy' 

y°/+  2j)!|  /”  *r"+A+2>  exp  (-r2)Mx(kBr).  (C.20) 


Knowing  that 


POO 

Q^(k,a)=  / drrN  exp  (—ar2)M\(kr) 
Jo 


(C.21) 


we  get 


A °o 


Q^(kA,kB,a)  = 


E 


( 


*iy 

2q  y 


«™^J!(2A  + l + 2j)!!^ 


g^+A+2j  / ! 


fci 
V^’ 


q(jv+a+i)/2  £ j\(2X  + l + 2j)\\Q'x  + + (C.22) 

We  only  need  to  calculate  Q^+A  and  Q^+x+2  the  others  are  obtained  by  recur- 
sion (type  1 integral  is  given  at  the  end  of  document) 


1 


Q?(k,  7)  = - 


k2  2n-5\  2 , 7i  -|-  4) (Z  n 3) 


7 L\47 


— + 


<?rJ+ 


q r 


. (C.23) 


The  upper  limit  to  the  utility  of  this  method  is  found  to  be  (A^  + fee)2 /2a  = 100 
when  approximately  70  terms  are  required  to  give  an  accuracy  of  10~13. 


Gaussian  Points  and  Weights  Method 

Writing  the  modified  spherical  Bessel  function  M\{z ) in  exponential  for  as 
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Mi(z)  = — [Ri{-z)exp(z)  - (-l)'i?/(z)exp(-2)] 


k=0 


k\(l  — k)\ 


(C.24) 

(C.25) 


For  large  z , Mi(z ) becomes  simply  the  first  term  in  the  Eq.  C.24.  Thus  when 
kA/ yfa  and  kB/y/a  are  large, 


Mt(z ) «=  —Ri(-z)  exp  (z);  R,(z)  « 1 « exp  (-z), 
ZZ  Zz 


(C.26) 


therefore 


roo 

Qyy(kA,kB,a)  — / drrN  exp  (—or2) 


exp  exp ( kBr ) 
2fc^r  2kBr 


r°° 

= — — — / drrN~2  exp  (— or2  + kAr  + fcBr).  (C.27) 
4A^b  Jo 

The  radial  integral  can  be  approximated  (after  a change  of  variables  r -A-  r/ y/a 


Q^{kA’kB'a) K L ir( 75)  exp{-r’+%r+%r)'ia2s) 

this  suggest  gaussian  integration.  Finding  the  maximum  by  differentiating  the  inte- 
grand 

N- 3 


2 , , kB 


(N  - 2)  ( -7=  ) "7=  exP  I -r'  + ~r  + — r 

\y/aj  y/a  \ y/a  y/a 


+ 


V® 


N- 2 


-2r  + = 0 (C.29) 

y/a  y/a ) 


2r  = r((^£±M!j+(Ar_2)  = 0 

1 (kA  + kB)2 


1 ( kA  + kB)  1 
rc  = 4 — 7S — =*=  5 


I 1/2 


a 


+ 2(iV  — 2) 


(C.30) 
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for  the  range  of  \{kA  + kB)/y/a  for  which  the  method  is  practical,  the  effect  of  the 
term  2 (N  — 2)  is  very  small.  Therefore,  keeping  rc  independent  of  N the  maximum 
is  approximated  as 


rc  = x- 


1 kA  + kB 


(C.31) 


2 ' 

A change  of  variables  t = r — rc  minimize  the  number  of  points  in  the  numerical 
integration: 


/OO 

dtf(t,  rc,  kA,  kB,  a)  exp  ( -t 2), 
rc 


(C.32) 


where 


f(t,  rc,  kA)  kB,  a)  = ~^=  M - 


x Mi 


kA 

y/a 

' kB 

.y/a 


( t + rc) 
(t  + rc) 


exp  [—2  ret  — r £].  (C.33) 


Extracting 


exp  (—4={t  + rc))  and  exp  {—{t  + rc)) 

y/a  y/a 

from  M\  and  Mx  in  Eq.  C.33  and  the  operating  exponents 


(C.34) 


„ o kA  , . kb 

%rct  -r%  + ~7={t  + rc)  + 
v“ 


hi 


n , 2 kA  + kB,  \ 

-2 rct  -rc  + + rc) 

y/a 

— >•  -2 rct  - r2c  + 2 rc(t  + rc)  — r2c  (C.35) 


then  / in  Eq.  C.33  reduces  to 


N 


f{t,rc,kA,kB,a)  = ~ J Mx 


_y/a 


{t  + rc) 


x 


M'x  + rc))  exp  (rc)’  (C-36) 
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where 


M'i  (z)  = ^[Ri(-z)-(-l)lRi(z)exp(-2z)].  (C.37) 

The  number  of  integration  points  required  for  a given  accuracy  decreases  with 
increasing  (A:^  + kB)2/2a. 


Range  of  ( kA  + kB)2/ 2a 

Number  of  points 

[102  ,103] 

20 

103,105 

10 

> 105 

5 

Finite  Sum 

This  method  can  be  applied  to  many,  but  not  all  of  the  needed  integrals.  A 
tabulated  integral  of  Bessel  function  in  the  present  notation  gives 


2 \ 7T1/2  /(H+A;d)\  ,,  ( kAkB 

Ql\(kA,  kB , a;)  = exp  I : — — 1 M; 


(C.38) 


(4a3/2)  V 4a  7 V 2 a 
By  differencing  this  expression  with  respect  to  kA  and  kB  various  numbers  of  times, 
using  the  Bessel  function  recursion  relation  and  using  the  binomial  coefficient  sum- 
mation relations  , a moderately  general  form  is  obtained: 

r1/2  /H+HX  X~Vti-V 


Q 


X+fi+2 — 2u 
A/i 


7T 


(4a3/2) 


exp 


4a 


t 


u 


EE 

t= 0 u=0 

(A  + fj,-v+  |)!(/i  + v - |)!(f  + u + v+\) 


A — iA  ( n — v 


(»+t  + m\  + u + Mv-\)\ 


X 


7 \ \ — U — t+U  / 7 \ /l  — I/+t  — U 

K4  \ / 


2a 


B 

2a 


M 


u-M+i' 


kAkB 

2a 


(C.39) 


To  avoid  underflows  and  overflows  a scaled  form  of  the  Bessel  function  is  used 

„A 


ma(z)  = (2A  + 1)!!  exp^)M^x)- 


(C.40) 
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We  then  obtain 


^\A+/i-f  2— 2*/ 
^A/i 


(2(A  + //  — v)  + 1)!! 


(2i/-  1)!! 


x exp 


/ (kA  + 

V (4a) 


A— i/  /i-i/ 

EE 

t=0  u=0 


A - ^ 


H — V 
, U 


(2{t  + v)  — l)!!(2(u  + u)  — 1)!!  (k\\u  ( kB V 
(2(/i  + t)  + 1)!!(2(A  + u)  + 1)!!  \2a ) V2aJ 


X MUu+u 


kAkB 
2 a 


1 

(2(t  + u + u)  - 1)!! 


where  we  have  used  M'  in  terms  of  Ms 


Mt(x ) = exp  (x)M[(x) 


x * 


(2/  + 1)!! 


exp  (x)Mi(x) 


_ (2^  + !)”  » w/  \ (2(A  + u + + 1)!!  ,,  v 


a;* 


(H^r 


The  integral  Q is  reduced  now  to 


(C.41) 

(C.42) 


^\A+/i+2— 2z/ 


A— i/ 


x exp 


/ (A^  + 

v 4 a 


kB 
2a 

\-v  H—V 

EE 

<=0  u=0 


(2(A  + n — u)  -T  1)!! 


A — u 


(2v  — 1)!! 


v 


u 


(2{t  + u)  - l)!!(2(u  + v)  - 1)!! 

(2(/i  + i)  + 1 ) ! ! (2 ( A + it)  + 1)!!  \kA ) 


x (2(t  + u + v)  + l)M't+u+v  ^ 2qB) 


(C.43) 


Type  1 Radial  Integral 


Q?(k,  a)  = y/nkl2  1 2a  (+2+  1 R{l,n)<f> 


(l  + n + l)  3 k2 
2 ’ + 2’4a 


(C.44) 
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here  R is  the  ratio  of  gamma  functions 


R(l,  n ) = 


r((/  + n + l)/2)_  (f  + n-1)!!, 
(2/  + 1)!! 


r(/  + §) 


2 2 +1  2 n + l even 


fi+"-iyi 

= V .2  o(*+1)  n + / odd 

y/n(2l  + 1)!!  ’ n + Zodd- 


(C.45) 

(C.46) 


In  Eq.  C.44  4>  is  the  degenerate  hypergeometric  function.  The  confluent  hyper- 
geometric series  for  (f)  is 


. . , , , a z a(a  + 1)  z2 

*M.e)  = l + -n  + ^L_i_  + 


(C.47) 
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